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NEW  COMPUTATIONAL  MODELS  IN  CONTINUUM  MECHANICS 


* 

O.M.  Belotserkovskii 

Direct  numerical  simulation  of  complex  gas  dynamics  problems  (compu- 
tational experiment)  is  performed  with  the  help  of  Euler,  Navier-Stokes 
and  Boltzmann  equations.  The  basic  principles  of  the  computational  experi- 
ment are  formulated  and  the  results  for  various  gas  dynamics  problems  of 
a complex  internal  structure  are  presented. 

The  problems  examined  include  the  transonic  regime  (super-critical 
flows  including  transition  through  sound  velocity),  flows  with  a jet 
"injected"  into  the  main  stream  and  diffraction  problems.  Body  wake  flows 
are  studied  at  various  Reynolds  numbers.  The  structure  of  a shock  wave  pro- 
vides an  example  of  rarefied  gas  flows  at  various  Mach  numbers. 

A set  of  control  tests  is  worked  out  for  the  estimation  of  calcula- 
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Preface 


Selected  results  of  this  report  and  its  companion,  "Investigation  of 
Transonic  Gas  Flows,"  were  communicated  in  seminar  lectures  given  by  Prof. 
Belotserkovskii  at  several  American  universities  ouring  a 4-week  visit  in 
November  - December,  1978.  His  lecture  on  "New  Computational  Models  in  Con- 
tinuum Mechanics,"  presented  December  1,  1978,  at  a seminar  of  the  Aerospace 
Engineering  Department  of  the  University  of  Maryland  at  College  Park,  at- 
tracted considerable  interest  from  a diverse  audience  of  scientists  and 
engineers.  In  response  to  numerous  requests  Prof.  Belotserkovskii  made 
available  the  manuscripts  of  two  reports  which  describe  in  some  detail  the 
computational  techniques  employed  in  the  numerical  solutions  of  the  problems 
surveyed  in  his  seminar  talks. 

This  report  is  a survey  of  computational  models  developed  by  the  author 
in  collaboration  with  his  colleagues  at  the  Computational  Center  of  the 
USSR  Academy  of  Sciences  over  the  past  decade.  Although  most  of  these  re- 
sults have  previously  appeared  in  various  Soviet  journals,  with  the  most 
recent  contribution,  the  last  chapter  of  the  present  report,  in  the  proceed- 
ings of  the  VI  International  Conference  of  Computational  Methods  of  Hydro- 
dynamics, Vol.  2,  Moscow,  1978,  pp.  37-47,  this  report  nevertheless  gives  a 
coherent  review  of  the  advances  in  computational  fluid  dynamics  at  one  of 
the  foremost  centers  of  the  Soviet  Union. 

The  hallmarks  of  their  numerical  techniques  are  that  they  are: 

Cl)  typically  differentially  and  globally  conservative,  and 

(2)  careful  attention  is  made  to  minimize  truncation  errors  while  at 
the  same  time  the  favorable  properties  of  the  computational 
schemes  are  exploited,  e.g.,  the  effective  viscosity  of  the 
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finite-difference  equations  (to  promote  calculational  stability) 
and  (physically  interpreted)  the  ability  to  compute  the  essential 
features  of  separated  regions  of  recirculating  flows  (e.g.,  wakes) 
wholey  within  the  framework  of  the  Euler  equations. 

The  underlying  theme  of  this  work  is  perhaps  best  described  by  the 
author.  "The  properties  that  a numerical  method  is  to  be  endowed  with, 
from  the  view-point  of  modem  developments,  are  so  diverse  that  they  are 
difficult  to  fully  implement  in  one  single  method.  In  view  of  this,  com- 
plexes of  numerical  methods  based  upon  a unified  approach  should  be  availa- 
ble. Finally,  it  is  desirable  to  consider  homogeneous  schemes  that  enable 
calculations  through  discontinuities  that  may  arise  in  the  evolution  of  the 
solution,  that  allow  for  explicitly  singling  out  some  (principal)  of  the 
features  and  that  adequately  resolve  their  boundary  conditions." 

The  numerical  methods  are  illustrated  by  a great  variety  of  computa- 
tional results  which  encompass  a wide  range  of  velocities,  from  subsonic 
through  transonic  as  well  as  up  to  hypersonic  wherein  complex  physical  pro- 
cesses (thermo- chemical  and  radiation)  strongly  affect  the  resulting  flow 
field  and  over  a wide  range  of  Reynold's  numbers  in  the  case  of  viscous 
compressible  gas  flows.  Wherever  possible  English  translations  of  the 
references  have  been  cited. 

The  careful  typing  of  the  edited  manuscript  by  Miss  Vicki  Brewer 
deserves  a special  note  of  appreciation.  Finally,  I am  pleased  to  acknow- 
ledge Contract  N0014-79-M-0022  from  the  U.S.  Office  of  Naval  Research  which 
made  possible  the  publication  and  distribution  of  this  report. 

W.  L.  Melnik,  Editor 

Professor,  Aerospace  Engineering  Dept. 

University  of  Maryland  at  College  Park 
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f 1 . Introduction 

% 

At  present,  specialists  of  applied  sciences  are  confronted  with 
various  kinds  of  practical  problems  whose  successful  and  accurate  solution, 
in  most  cases,  may  be  obtained  only  by  numerical  methods  with  the  aid  of 
computers . However,  this  does  not  mean  that  analytical  methods  which 
enable  us  to  find  the  solution  in  "closed"  form  will  not  be  developed. 
Nevertheless,  it  is  absolutely  clear  that  the  range  of  problems  permitting 
such  an  approach  to  their  solution  is  rather  narrow,  therefore,  the  de- 
velopment of  general  numerical  algorithms  for  the  investigation  of  prob- 
lems of  mathematical  physics  (gas  dynamics,  theory  of  elasticity,  etc.) 
is  especially  important. 

a)  Difficulties  of  Carrying  Out  the  Experiment. 

At  hypersonic  flight  velocities,  the  resulting  high  temperatures  give 
rise  to  effects  of  dissociation  and  ionization  in  the  flow  and,  in  a num- 
ber  of  cases,  even  to  "luminescence"  of  the  gas.  In  these  cases  it  is 
enormously  difficult  to  simulate  the  experiment  in  the  laboratory,  since 
it  is  not  sufficient  to  satisfy  the  classical  criteria  of  similarity, 
i.e.,  the  equality  of  the  Mach  and  Reynolds  numbers . The  equality  of  abso- 
lute pressures  and  absolute  temperatures  is  also  required,  which  is  only 
possible  if  the  sizes  of  the  model  and  prototype  are  identical.  These 
requirements  pose  numerous  technical  difficulties  and  not  the  least  of 
which  concerns  the  high  cost  of  the  experiment. 

However,  the  importance  of  the  experiment  must  not  be  underestimated 
for  it  is  always  the  basis  of  measure  for  confirming  (or  rejecting)  the 
calculation  scheme  and  numerical  solution. 
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b)  Complexity  of  the  Equations  Considered. 

The  prominance  of  numerical  methods  in  mechanics  of  continua  is  partly 
explained  by  the  fact  that  the  governing  equations  of  aerodynamics  and 
gas  dynamics,  and  of  theory  of  elasticity  represent  the  most  complicated 
(as  compared  to  other  branches  of  mathematics)  system  of  partial  differen- 
tial equations.  In  the  general  case,  this  is  a nonlinear  system  of  mixed 
type  (wherein  the  surface  across  which  the  equations  change  their  type  is 
unspecified)  and  with  "movable  boundaries",  i.e.,  the  boundary  conditions 
are  given  on  surfaces  or  lines  which,  in  turn,  are  determined  by  the 
calculations.  Moreover,  the  range  of  the  unknown  functions  is  so  wide 
that  ordinary  methods  of  analytical  investigation  (linearization  of 
equations,  series  expansion,  separation  of  a small  parameter,  etc.)  are 
not  appropriate  for  the  development  of  the  complete  solution  of  the  pro- 
blem in  the  general  case. 

It  should  be  noted,  that  in  solving  complicated  problems  on  electronic 
computers  that  the  preliminary  analytical  investigation  of  a problem  may 
provide  great  insight  and  sometimes  this  investigation  is  simply  decisive 
for  the  successful  realization  of  the  numerical  algorithm. 

Let  us  dwell  on  one  more  peculiarity  of  algorithms  used  for  solving 
concrete  problems  of  mechanics  of  continua.  Currently,  numerical  methods 
have  found  a wide  practical  application  in  design  offices  and  research 
institutes.  Substantial  progress  in  the  exploration  of  the  cosmos,  in 
the  optimum  control  of  vehicles,  in  the  choice  of  rational  configurations 
of  vehicles  and  etc.,  are,  to  a considerable  extent,  due  to  scientific 
information  obtained  from  serial  calculations . The  volune  of  information 
obtained  by  means  of  the  calculation  is  far  more  complete  and  substantially 


2 


cheaper  than  the  corresponding  experimental  investigations  if  the  problem 
is  correctly  formulated,  well  simulated  and  algorithmically  rational.  How- 
ever, a wide  application  of  the  numerical  methods  for  practical  purposes 
requires  sufficient  simplicity  and  reliability.  Thus,  on  the  one  hand,  one 
has  to  deal  with  rather  complicated  mathematical  problems,  while  on  the 
other  hand,  it  is  necessary  to  develop  rather  simple  and  reliable  numeri- 
cal methods  permitting  us  to  carry  out  serial  calculations  at  project 
institutes  and  design  offices. 

Note  that  for  most  problems  in  gas  dynamics,  not  only  have  no  mathe- 
matical theorems  of  existence  and  uniqueness  been  proved,  but  very  often 
there  is  no  confidence  that  such  theorems  can  even  be  derived.  As  a rule, 
the  veiy  mathematical  formulation  of  the  problem  is  not  strictly  given 
and  only  the  physical  treatment  is  presented,  which  is  far  from  being  one 
and  the  same  thing.  The  mathematical  difficulties  of  the  investigation  of 
such  types  of  problems  are  related  to  the  nonlinearity  of  the  equations, 
as  well  as  to  the  number  of  independent  variables. 

The  state  of  affairs  with  the  methods  of  solution  are  no  better.  So 
far,  investigations  related  to  the  possibility  of  realization  of  the  algo- 
rithm, its  convergence  to  the  unknown  solution,  and  its  stability  have 
rigorously  been  preformed  only  for  linear  systems,  and,  in  a number  of 
cases,  only  for  equations  with  constant  coefficients.  When  confronted 
with  the  necessity  of  solving  a problem,  the  mathematician  often  has  to 
use  the  known  algorithms  and  to  develop  new  methods  without  a rigorous 
mathematical  basis  for  their  applicability. 

In  science,  as  well  as  in  mathematics,  one  can  find  many  examples 
when  new  ideas  and  concepts  were  successfully  used  without 
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a rigorous  basis  which  only  appeared  later.  Of  course,  this  does  not 
suggest,  that  when  developing  new  numerical  algorithms,  one  may  slight 
the  accurate  formulation  of  the  problem  or  its  physical  meaning.  This 
oversight  inevitably  leads  to  numerous  mistakes,  consequently  a waste  of 
time  and,  moreover,  the  experience  without  being  theoretically  interpreted 
does  not  give  the  foundation  for  further  development  of  the  method. 

We  want  to  draw  your  attention  to  this  rather  clear  question  only 
because  there  is  still  prevalent  an  opinion,  that  the  main  thing  is  to 
write  down  differential  equations  and  all  the  rest  reduces  to  a trivial 
substitution  of  derivatives  by  differences  and  to  programming  on  which 
too  much  importance  is  sometimes  attached.  In  this  connection,  it  is 
reasonable  to  formulate  the  main  stages  of  the  numerical  solution  to  a 
physical  problem  with  the  aid  of  computers  in  the  following  way: 

1)  the  construction  of  a physical  model  and  the  mathematical  state- 
ment of  the  problem; 

2)  the  development  of  a numerical  algorithm  and  its  theoretical 
interpretat ion ; 

3)  programming  (manual  or  automatic)  and  the  formal  adjustment  of 
the  program; 

4)  the  methodical  adjustment  of  the  algorithm,  i.e.,  the  test  of 
its  operation  by  concrete  problems;  the  elimination  of  drawbacks  uncovered 
and  the  experimental  investigation  of  the  algorithm; 

5)  accumulation  of  experience,  the  estimation  of  effectiveness  and 
the  range  of  applicability  of  the  algorithm  from  serial  calculations. 

At  all  stages,  the  mathematical  theory,  the  physical  and  computational 
experiment  are  used  jointly  and  consistently.  Their  application  may  be 
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illustrated  by  solving  concrete  problems,  which  will  be  described  below. 
Therefore,  we  shall  make  only  some  common  observations  by  way  of  intro- 
duction. 

The  main  principle  of  using  mathematical  results  is  that  the  condi- 
tions providing  the  solution  of  a problem  for  special  simplified  cases 
must  be  fulfilled  as  well  for  more  general  cases.  Parallel  to  this,  con- 
sideration of  the  physical  phenomenon  provides  a qualitative  picture 
against  which  the  statement  of  the  problem  is  checked  and  defined  more 
exactly.  Ultimately,  the  final  experimental  test  allows  us  to  access  the 
correctness  of  the  assumptions  made  and  the  validity  of  the  algorithm  and 
resulting  solution.  It  should  be  noted,  that  the  estimate  of  accuracy  of 
the  numerical  solution  must  be  done  purely  mathematically,  without  using 
the  results  of  the  physical  experiment.  The  latter  may  be  used  for  quali- 
tative comparison  while  the  quantitative  comparison  between  the  calcula- 
tion and  experiment  provides  information  on  how  closely  the  physical  model 
used  approaches  the  natural  environment. 

1.1  Numerical  Methods  of  Solution  of  Equations  of  Gasdynamics. 

Many  important  problems  of  the  exact  sciences  involve  the  solution  of 
a system  of  non-linear  partial  differential  equations.  Oftentimes  it  includes 
many  problems  with  discontinuous  solutions  (e.g.,  gas  dynamics). 

The  construction  of  reasonably  accurate  solutions  of  the  exact  equa- 
tions of  gasdynamics  in  the  general  case  has  become  possible  only  with  the 
aid  of  mmerical  methods,  exploiting  the  advantages  of  high-speed  elec- 
tronic digital  computing  machines.  Technological  requirements  have  called 
for  an  intensive  development  of  mmerical  methods  and  their  application  to 
the  solution  of  a wide  variety  of  gas  dynamics  problems.  Scientists  and 
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research  engineers  in  the  area  of  gas  dynamics  have  contributed  signifi- 
cantly to  the  development  of  modem  numerical  methods  of  solving  systems 
of  non-linear  partial  differential  equations. 

There  exist  four  universal  numerical  methods  which  are  applicable  to 
the  solution  of  non-linear  partial  differential  equations  of  gas  dynamics 
problems . 

J 

I.  Method  of  Finite  Differences.  This  method  is  the  most  highly 
developed  of  the  four  at  the  present  time  and  is  widely  applied  to  the 
solution  of  both  linear  and  non-linear  equations  of  the  hyperbolic,  ellip- 
tic and  parabolic  types.  The  region  of  integration  is  subdivided  into  a 
network  of  computational  cells  by  a generally  fixed  orthogonal  mesh. 
Derivatives  of  functions  in  the  various  directions  are  replaced  by  finite 
differences  of  one  form  or  another;  usually,  a so-called  implicit  differ- 
ence scheme  is  applied  to  the  integration  of  the  equations.  This  results 
in  the  solution,  at  each  step  of  the  procedure,  of  a system  of  linear 
algebraic  equations  involving  perhaps  several  hundred  unknowns. 

Finite  difference  schemes  are  often  used  for  solving  steady  and  un- 
steady gas  dynamics  equations.  Lagrangian  and  Eulerian  approaches  are 
widely  used  here.  In  the  first  case,  where  the  coordinate  network  is  re- 
lated to  the  fluid  particles  the  structure  of  the  flow  is  better  defined 
and  one  succeeds  in  constructing  rather  accurate  numerical  schemes  for 
flows  with  comparatively  small  relative  displacements.  In  the  second 
case,  when  the  calculational  network  is  fixed  over  space,  the  schemes  are 
used  for  constructing  flows  with  large  deformation.  In  recent  time,  the 
approaches  mentioned  here  have  also  found  a wide  application  to  the  cal- 
culation of  steady  flows. 
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II.  Method  of  Integral  Relations.  In  this  method,  which  is  a 
generalization  of  the  well-known  method  of  straight  lines,  the  region  of 
integration  is  subdivided  into  strips  by  a series  of  curves,  the  shape  of 
which  is  determined  by  the  form  of  boundaries  of  the  region.  The  system 
of  partial  differential  equations  written  in  divergence  form  is  integrated 
across  these  strips,  the  functions  occuring  in  the  integrands  being  re- 
placed by  known  interpolation  functions.  Hie  resulting  approximate  system 
of  ordinary  differential  equations  is  integrated  numerically.  The  method 
of  integral  relations,  like  the  method  of  finite  differences,  is  applica- 
ble to  equations  of  various  types. 

III,  Method  of  Characteristics.  This  method  is  usually  only  applied 
to  the  solution  of  equations  of  hyperbolic  type.  The  solution,  in  this 
case,  is  computed  with  the  aid  of  a grid  of  characteristic  lines,  which 

is  constructed  in  the  course  of  the  computation.  Actually,  the  method  of 
characteristics  is  a difference  method  of  integrating  systems  of  hyperbolic 
equations  on  the  characteristic  calculational  network  and  is  mainly  used 
for  detailed  description  of  flows.  Its  distinguishing  feature  as  compared 
to  other  difference  mentods  is  the  minimal  utilization  of  interpolation 
operators  and  associated  with  it  the  region  of  influence  of  the  difference 
scheme  closely  approximating  the  region  of  influence  of  the  system  of 
differential  equations.  The  smoothing  of  the  profiles,  inherent  in 
difference  schemes  with  fixed  network,  is  minimized  since  the  calculational 
network  used  in  the  method  cf  characteristics  is  constructed  exactly 
with  the  region  of  influence  of  the  system  taken  into  account. 

Irregularity  (nonconservativeness)  of  the  calculational  network 
should  be  noted  as  a drawback  of  the  method  of  characteristics.  It  is 
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possible  to  develop  a technique,  based  on  this  method,  in  which  the  calcu- 
lations are  carried  out  in  layers  bounded  by  fixed  lines.  The  method  of 
characteristics  permits  one  to  accurately  determine  the  point  of  origin 
of  secondary  shock  waves  within  the  field  of  flow  as  the  result  of  inter- 
section of  characteristics  of  oijie  family.  On  the  other  hand,  if  a large 
number  of  such  shock  waves  occur,  difficulties  are  encountered  in  their 
calculation.  Accordingly,  the  method  of  characteristics  is  most  expedient- 
ly applied  to  hyperbolic  problems  in  which  the  number  of  discontinuities 
is  small  (for  example,  problems  concerning  steady  supersonic  gas  flow). 

IV.  "Particle-in-Cell"  Method  (PIC).  In  certain  respects,  the  PIC 
method  incorporates  the  advantages  of  both  the  Lagrangian  and  Eulerian 
approaches.  The  range  of  solution  here  is  separated  by  the  fixed  (Eulerian) 
calculation  network;  however,  the  continuous  medium  is  represented  by  a 
discrete  model,  i.e.,  the  population  of  particles  of  fixed  mass  (Lagran- 
gian network  of  particles)  which  move  across  the  Eulerian  network  of  cells 
is  considered.  The  particles  are  used  to  determine  parameters  of  the 
fluid  itself  (mass,  energy,  velocity),  whereas  the  Eulerian  network  is 
employed  for  determining  parameters  of  the  field  (pressure,  density, 
temperature) . 

The  PIC  method  allows  us  to  investigate  complex  phenomena  of  multi- 
component media  in  dynamics,  because  particles  carefully  monitor  free 
surfaces,  lines  of  separation  of  the  media,  etc.  Due  to  discrete  repre- 
sentation of  a continuous  medium  (the  finite  number  of  particles  in  a cell) 
calculational  instability  (fluctuations)  often  occurs.  However  the  cal- 
culation of  rarefied  regions  is  also  difficult.  Limitations  in  capacity 
of  modem  computers  do  not  permit  a significant  . tcrease  in  the  number  of 
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particles. 

For  problems  in  gas  dynamics  involving  a uniform  mediun,  it  seems 
more  reasonable  to  employ  the  concept  of  continuity  considering  the  mass 
flow  across  the  boundaries  of  Eulerian  cells  instead  of  "particles". 

Only  numerical  methods  using  high  speed  computers  and  careful  experi- 
ments allow  us  to  obtain  the  complete  solution  to  complex  gas  dynamics 
problems  and  to  determine  the  necessary  flow  characteristics.  Thus,  the 
elaboration  of  numerical  schemes,  the  calculation  of  different  gas  dyna- 
mics problems,  as  well  as  the  study  of  analytical  properties  of  the  solu- 
tions and  their  asymptotic  behavior  are  of  significant  interest  at  present. 

1.2  Aerodynamic  Characteristics  of  High-Speed  Vehicles. 

In  this  paper  a short  review  of  the  numerical  methods  used  for  the 
determination  of  the  aerodynamic  characteristics  of  high-speed  vehicles 
with  transonic  and  supersonic  velocities  will  first  be  given.  The  numeri- 
cal schemes  were  developed  under  our  supersivion  and  in  collaboration  witli 
the  Moscow  Physical  Technical  Institute  and  the  Computing  Center  of  the 
Academy  of  Sciences  of  the  USSR.  We  shall  discuss  the  problems  of  the 
development  and  use  of  the  numerical  algorithms  for  carrying  out  serial 
calculations  in  solving  modem  engineering  problems  arising  in  practice. 

I.  Steady -State  Schemes.  In  determining  the  steady  aerodynamic 
characteristics  of  bodies  (especially  when  electronic  computers  of  average 
capacity  were  employed)  we  made  wide  use  of  the  following  methods  for  sol- 
ving the  steady  gas  dynamics  equations:  the  method  of  integral  relations 
(m.i.r.),  the  method  of  characteristics  (m.ch.)  and  some  finite  difference 
schemes  (e.g.,  schemes  with  "artificial  viscosity").  We  wish  to  especially 
consider  problems  in  which  different  discontinuities  and  singularities  are 
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given  beforehand,  together  with  some  associated  boundary  conditions;  the 
solutions  being  carried  out  in  regions  where  functions  vary  continuously. 

Three  different  schemes  of  the  method  of  integral  relations  have  been 
developed  for  the  determination  of  flow  in  the  nose  region  of  a blunt  body, 
namely,  one  that  employs  an  interpolation  of  various  functions  across  the 
shock  layer  (Scheme  I),  along  it  (Scheme  II)  or  in  both  directions  (Scheme 
III) . As  a result  the  boundary  value  problem  is  solved  for  an  approxi- 
mate system  of  ordinary  differential  equations  (Schemes  I and  II)  or 
algebraic  equations  (Scheme  III) . To  solve  the  three  dimensional  problem, 
certain  additional  trigonometric  approximations  in  the  circumferential 
coordinate  were  introduced.  The  various  schemes  of  the  method  of  integral 
relations  have  found  a wide  variety  of  applications  [1-6]. 

The  main  advantage  of  these  schemes  is  that,  by  means  of  different 
transformations,  one  succeeds  eventually  in  approximating  functions  (or 
groups  of  functions)  with  comparatively  weak  variations.  Consequently 
reliable  results  with  a high  degree  of  accuracy  can  be  obtained  with  a 
comparatively  small  number  of  interpolation  nodes  (usually  3-4  are  used) . 

The  choice  of  the  independent  variables,  the  form  of  the  initial 
system  of  equations  of  motion  (i.e.,  the  introduction  of  the  integrals 
into  the  initial  system  and  the  use  of  the  divergence  form  of  the  laws  of 
conservation) , the  use  of  conservation  schemes , the  approximation  of  the 
integrals,  etc.,  are  all  of  great  importance  in  writing  an  efficient  numeri- 
cal algorithm  using  m.i.r.  . 

The  main  difficulty  in  carrying  out  the  schemes  of  m.i.r.  is  the 
solution  of  many  parameter  boundary  value  problems  for  the  approximating 
system  of  equations.  This  is  overcome  by  means  of  appropriate  iteration 
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schemes.  Moreover,  these  schemes  have  been  used  in  transonic  regions 
mainly  for  bodies  of  a comparatively  simple  form,  while  when  dealing  with 
a supersonic  zone  one  has  to  adopt  another  algorithm. 

In  calculating  supersonic  flow  the  two-  and  three-dimensional  schemes 
of  the  method  of  characteristics  by  P.I.  Qiushkin,  K.M.  Magomedov,  and 
their  co-workers  were  used  [7,8].  With  the  governing  equations  expressed 
in  characteristic  variables,  one  requires  approximation  of  ordinary  deri- 
vatives, only.  Using  a fixed  computational  network,  a system  of  linear 
finite  difference  equations  is  obtained  with  its  attendant  advantages. 

With  the  help  of  the  methods  cited  above,  a large  number  of  gas 
dynamics  problems  have  been  solved,  namely,  ideal  gas  flows  with  chemical 
reactions  and  radiation,  transonic  and  three-dimensional  flows,  as  well 
as  viscous  flows.  In  most  cases  reliable  results  were  obtained  which 
were  in  excellent  agreement  with  experiment  [6].  However,  these  approaches 
to  the  solution  of  the  steady- state  equation  may  be  successfully  used  only 
for  problems  in  which  there  are  no  singularities,  discontinuities,  inter- 
sections, and  interactions.  The  application  of  these  approaches  is  difficult 
for  bodies  of  complex  form  with  a large  number  of  discontinuities.  Besides, 
a single  algorithm  for  the  calculation  of  different  types  of  flow  is  pre- 
ferable. 

II.  Unsteady-State  Schemes.  The  next  step  in  the  evolution  of  nuneri- 
cal  methods,  which  was  motivated  by  urgent  practical  needs  and  aided  by 
the  availability  of  electronic  computers , was  the  development  of  nonsteady 
schemes  to  calculate  the  asymptotic  solution  of  steady- state  aerodynamic 
problems.  In  approximating  the  nonsteady  equations,  the  general  principles 
and  ideas  of  the  m.i.r.  and  m.ch.  were  applied  with  respect  to  space 
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variables.  The  divergence  or  characteristic  forms  of  the  initial  equa- 
tions were  used  and  the  same  calculational  networks  were  employed. 

In  this  way  the  nonsteady  Schemes  II  and  III  of  the  method  of  inte- 
gral relations  and  the  network-characteristic  method  were  developed  [9,10]. 
In  this  way  rather  complicated  types  of  flow  could  be  treated  with  a sin- 
gle algorithm.  It  is  natural  that  the  problems  of  computational  stability 
and  the  attainment  of  steady- state  solutions  should  become  crucial.  They 
require  some  specific  technique  such  as  the  introduction  of  artificial 
viscosity  into  the  initial  system,  and  of  dissipation  terms  into  the 
difference  equations.  In  a number  of  cases  the  accuracy  of  the  results 
obtained  is  somewhat  less  than  in  the  steady- state  methods,  but  these 
approaches  enabled  us  to  consider  new  classes  of  problems;  for  example, 
the  determination  of  the  aerodynamic  characteristics  of  specific  configura- 
tions involving  three-dimensional  flow,  the  calculations  of  viscous  tran- 
sonic flows,  and  others  [10]. 

III.  "Large  Particles"  Method.  Finally,  in  the  third  stage  of 
development  it  seemed  reasonable  and  advantageous  to  introduce  the  ele- 
ments of  the  Harlow  "particle -in- cell"  method  [11-13]  into  the  algorithms. 
At  first  only  the  equation  of  continuity  is  represented  as  the  mass  flow 
across  the  Euler  cell,  using  the  simplest  finite  difference  or  integral 
approximation  along  the  coordinates. 

Thus  the  modified  method  of  "large  particles"  [14-15]  came  into  exis- 
tence, which  (again  by  means  of  the  stability  process)  allowed  us  to  con- 
sider from  one  point  of  view  such  a complicated  task  as,  for  example,  the 
supsonic,  transonic,  and  supersonic  flow  past  a flat-nosed  body  in  two 
dimensions  or  with  axial  symmetry.  This  method  is  likewise  used  in  cal- 
culating viscous  flows  which  would  even  permit  the  study  of  separated  flows 
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from  solutions  of  the  Navier-Stokes  equations  [16-18]. 

. The  main  principle  of  splitting  the  evolution  of  a physical  process 

by  a time  step  is  as  follows. 

The  medium  simulated  may  be  replaced  by  a system  of  N particles 
(fluid  particles  for  a continuous  medium  and  molecules  for  a discrete  one) 
which  at  the  initial  instant  of  time  are  distributed  in  cells  of  the 
Eulerian  mesh  in  a coordinate  space  in  accordance  with  the  initial  data. 

The  evolution  of  such  a system  in  time  At  may  be  split  into  two  stages: 
change  of  the  internal  state  of  sub-systems  in  ceils  which  are  assumed  to 
be  "frozen"  or  stable  ("Eulerian"  stage  for  a continuous  medium  and  colli- 
sion relaxation  for  a discrete  one)  and  subsequent  displacement  of  all  the 
particles  proportional  to  their  velocity  and  At  without  changing  the  inter- 
nal state  ("Lagrangian"  stage  for  a continuous  medium  and  free  motion  of 
molecules  for  a discrete  one) . The  stationary  distribution  of  all  the 
medium  parameters  is  calculated  after  the  process  attains  steady  state. 

It  should  be  stressed  that  the  development  of  the  numerical  schemes 
mentioned  above  has  been  paced  by  the  improvement  and  extension  of  the 
ways  of  solving  the  boundary  valuie  problems  for  the  corresponding  approxi- 
mating equations;  by  the  consideration  of  a new,  wider  class  of  problems; 
by  the  development  and  improvement  of  electronic  computers,  machine 
languages,  input  and  output  arrangements,  etc. 

1.3  Computational  Experiments 

In  recent  years  the  introduction  of  big  computers  has  aroused  a con- 
siderably greater  interest  in  various  numerical  methods  and  algorithms 
whose  realization  borders  on  carrying  out  computational  experiments.  The 
• need  in  such  an  approach  for  the  solution  of  problems  of  mathematical 
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physics  is  prompted  by  ever  growing  practical  demands;  in  addition  it  is 
connected  with  an  attempt  of  constructing  more  rational  general  theore- 
tical models  for  the  investigation  of  complex  physical  phenomena. 

Let  us  outline  the  principle  steps  of  computational  experiments.  At 
first,  one  chooses  a mathematical  model  of  a physical  object  based  on  its 

analytical  study.  Then  one  works  out  a tool  for  the  investigation  of  the 

phenomenon  in  question,  namely  a difference  scheme  which  permits  the  ex- 
periment itself  to  be  carried  out,  i.e.,  the  computational  process.  The 
next  step  comprises  a detailed  analysis  of  the  results,  leading  to  improve- 
ments and  corrections  of  the  mathematical  model.  This  feedback  procedure 
leads  to  perfections  and  modifications  in  the  methodology  of  numerical 
experiments . 

A close  analogy  to  physical  experiments  comprising  similar  steps  is 

evident;  an  analysis  of  the  phenomenon  under  study;  development  of  an 

experimental  scheme;  modification  of  design  elements  of  the  experimental 
installation;  and  measurements  and  their  analysis. 

In  recent  years  the  Computing  Center  of  the  USSR  Academy  of  Sciences 
carried  out  a number  of  experiments  associated  with  studies  of  complex 
gasdynamic  flows  using  the  non-stationary  method  of  "large  particles" 
[14,15].  Characteristic  features  of  flows  past  bodies  of  different  shapes 
were  studied  over  a wide  range  of  velocities,  from  subsonic,  through  tran- 
sonic, up  to  hypersonic.  In  this  paper  results  of  a number  of  such  experi- 
ments are  presented  without  delving  into  the  details  of  the  computations 
[10,14,15,19,20]. 

It  also  seemed  promising  to  apply  the  main  principles  of  the  approach 
in  question  for  the  simulation  of  rarefied  gas  flows.  The  application  of 
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a statistical  variant  of  such  an  approach  for  the  solution  of  the  Boltz-  \ 

mann  equation  is  studied  in  [21-23].  Since  the  complete  details  of  the 

techniques  are  given  in  the  cited  references,  this  paper  will  only  be  -\ 

I 

concerned  with  the  characteristic  features  of  each  approach.  \ 
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2.  "Large  Particles"  Method  for  the  Study  of  Complex  Gas  Flows. 

For  numerical  models  constructed  by  Yu.M.  Davidov  [14,15]  on  the  basis 
of  Ruler ian  equations,  the  mass  of  a whole  fluid  (Eulerian)  cell,  i.e., 

"a  large  particle"  (from  which  originates  the  name  of  the  method)  is  con- 
sidered instead  of  the  ensemble  of  particles  in  cells.  Furthermore,  non- 
-stationary  (and  continuous)  flows  of  these  "large  particles"  across  the 
Eulerian  mesh  are  studied  by  means  of  finite-difference  or  integral  re- 
presentations of  conservation  laws. 

This  method  utilizes  the  conservation  laws  given  in  the  form  of  ba- 
lance equations  for  a cell  of  finite  dimensions  (which  is  the  usual  proce- 
dure in  deriving  the  gas  dynamics  equations  but  stops  short  of  passing  in 
the  limit  to  a point.  As  a result,  we  obtain  divergence -form  conservative 
and  dissipative-steady  numerical  schemes  that  allow  us  to  study  a wide 
class  of  complex  gas  dynamics  problems  (transonic  flows,  turbulent  flows 
in  the  wake  of  a body,  diffraction  problems,  etc.)  [10,14,15,20]. 

2.1  Description  of  the  Method 

Consider  the  motion  of  an  ideal  compressible  gas.  Our  starting-point 
is  provided  by  the  Euler  differential  equations  in  divergence  form  (the 
equations  of  continuity,  momentum  and  energy) 

f£+  V • (p$)  = 0, 

v.fprfwfc  .0.  m 

v . c»v\b  * f£  - 0, 

+ V • (pEV)  + V • (pV)  = 0. 

It  was  shown  in  [14]  that,  in  the  "large  particle"  method,  the  set  of  equa- 
tions of  gas  dynamics,  written  as  laws  of  conservation  in  integral  form, 
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may  be  used  instead  of  Cl).  The  important  point  is  that  the  difference 
scheme  approximating  the  initial  set  of  equations  should  be  homogeneous, 
so  that  "through"  computations  may  be  performed  without  isolating  singu- 
larities. 

Equations  (1)  are  completed  by  the  equation  of  state 

p = p(p,E,$)  = 0.  (2) 


The  various  stages  of  the  computational  cycle  will  be  considered 
separately.  Let  us  briefly  describe  the  main  principles  of  the  "large 
particle"  method.  The  region  of  integration  is  covered  by  a fixed  (over 
space)  Euler  mesh  composed  of  rectangular  cells  with  sides  Ax,  Ay  (or  az, 
Ar  in  a cylindrical  coordinate  system). 

In  the  first  ("Eulerian")  stage  of  calculations  only  those  quantities 
change  which  are  related  to  a cell  as  a whole,  and  the  fluid  is  supposed 
to  be  momentarily  decelerated.  Hence,  the  convective  terms  of  the  form 
div($pV)  where  4>  = (l,u,v,E),  corresponding  to  displacement  effects,  are 
omitted  in  equation  (1) . Then  it  follows  from  the  equation  of  continuity, 
in  particular,  that  the  density  field  will  be  "frozen"  and  the  initial 
Systran  of  equations  will  be  of  the  form 
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Here  we  have  used  both  the  simplest  finite -difference  approximations 
and,  to  improve  the  calculation  stability,  the  schemes  of  the  method  of 
integral  relations,  in  which  "sweeping-through"  approximations  of  the  in- 
tegrands with  respect  to  rays  (N  = 3,4,5)  are  used. 

In  the  second  ("Lagrangian")  stage  we  find  mass  flows  across  the 
cell  boundaries  at  tn+^  = tn  + At.  At  this  stage  we  assume  total  mass  to 
be  transferred  only  by  a velocity  component  normal  to  the  boundary.  Thus, 
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for  instance 


^,3  ' <uU,j>ly  11 
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The  angle  brackets  ^ ^ denote  the  value  of  p and  u on  the  cell  boundary. 

The  choice  of  these  values  is  extremely  important  since  they  substantially 
influence  the  stability  and  accuracy  of  the  computation.  The  various  pos- 
sible ways  of  writing  down  AM11  are  characterized  by  consideration  of  the 
flow  direction. 

First  and  second  order  accurate  representations  of  AM11  are  considered. 
These  are  based  on  central  differences,  without  account  being  taken  of  the 
flow  direction,  as  well  as  by  means  of  the  discrete  model  of  a continuous 
medium  comprising  a combination  of  particles  of  a fixed  mass  in  a cell 
[14,15]. 

Lastly,  in  the  third  ("Final")  stage  we  estimate  the  final  fields  of 
the  Euler  flow  parameters  at  the  instant  of  time  tn+1  (all  the  errors  in 
the  solution  of  equations  are  "removed") . As  was  pointed  out,  the  equa- 
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tions  at  this  stage  are  laws  of  conservation  of  mass  M,  momentum  P and 
total  energy  written  down  for  a particular  cell  in  the  difference  form 

Fn+1  = F11  + I A V*dTy  , where  F = (M,?,E)  (5) 

According  to  these  equations,  inside  the  flow  field  there  are  no 
sources  or  sinks  of  M,  ? and  E and  their  variation  in  time  At  is  caused 
by  interaction  at  the  external  boundary  of  the  flow  region. 

2.2  Boundary  Conditions. 

To  retain  the  unified  nature  of  the  computations  and  avoid  special 
expressions  for  the  boundary  cells,  layers  of  fictitious  cells  are  intro- 
duced along  all  the  boundaries,  which  are  assigned  parameters  from  the 
neighboring  flow  cells.  The  mrnber  of  such  layers  depends  on  the  order 
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of  the  difference  scheme  (one  layer  for  the  scheme  of  first-order  accuracy , etc.). 
Two  kinds  of  boundary  then  have  to  be  distinguished:  the  rigid  boundary 
(or  axis  of  symmetry)  and  the  "open"  boundary  of  the  computational  region. 


In  the  first  case,  the  velocity  component  normal  to  the  boundary 
changes  sign,  i.e.,  non-penetration  condition  along  rigid  walls,  while  the 
remaining  flow  parameters  are  taken  unchanged.  However  another  type  of 
boundary  condition  is  possible,  namely,  walls  without  slip  (condition  of 
sticking) . In  this  case  both  velocity  components  change  sign  and  the 
entire  velocity  vector  vanishes  on  the  wall. 

Fluid  can  flow  across  "open"  boundaries  of  the  region,  and  some  con- 
ditions on  the  continuity  of  the  movement  are  required  in  this  case. 
Gonsider  the  fluid  to  be  flowing  into  the  rectangular  mesh  from  the  left; 
then  the  parameters  of  the  entering  flow  will  be  specified  here.  On  the 
remaining  "open"  boundaries  of  the  region  we  extrapolate  the  parameters 
of  the  flow  "from  within",  i.e.,  transfer  to  the  fictitious  layer  the 
parameter  values  from  the  layer  nearest  to  the  boundary  (zero  order  extra- 
polation). A more  complicated  statement  of  the  conditions  is  possible, 
or  more  accurate  extrapolation  (say  linear  or  quadratic). 

It  is  natural  that  the  outer  boundary  of  the  region  should  be  fairly 
remote  from  the  source  of  disturbance,  in  which  case  methods  of  flow  extra- 
polation "outwards"  are  possible.  This  topic  will  be  discussed  in  more 
specific  terms  below.  It  will  merely  be  mentioned  here  that  the  basic 
principle  underlying  the  statement  of  the  conditions  is  that  no  substantial 
disturbances  should  penetrate  through  the  "open"  boundaries  of  the  region 
into  the  computational  region. 

2.3  Viscosity  Effects. 

It  has  already  been  remarked  that  this  approach  employs  homogeneous 
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difference  schemes,  whereby  computation  by  a unified  algorithm  is  possible 
both  through  smooth  flow  regions  and  through  discontinuities.  This  is 

achieved  by  using  finite-difference  schemes  with  a viscosity  approximation. 
Let  us  dwell  briefly  on  this  topic. 

While  the  equations  of  gas  dynamics  for  a non-viscous  gas  were  taken 
as  the  governing  equations,  viscosity  effects  are  in  fact  inherent  in  our 
difference  scheme.  They  are  produced,  firstly,  by  the  introduction  into 
the  scheme  of  an  explicit  term  with  artificial  viscosity  ('Viscosity  pres- 
sure") and  secondly,  by  the  presence  of  an  essentially  schematic  viscosity, 
dependent  on  the  structure  of  the  finite-difference  equations. 

The  form  of  the  approximation  viscosity  and  estimates  for  the  stabili- 
ty of  the  scheme  can  be  obtained  by  writing  as  Taylor  series  the  difference 
operators  appearing  in  the  equations  in  all  three  stages.  The  terms  of 
zero  (lowest)  order  should  then  represent  the  initial  differential  equa- 
tions, while  the  structure  of  the  approximation  viscosity  can  be  determined 
by  retaining  higher  order  terms  in  the  expansions  ("expansion  errors") . 

The  resulting  differential  equations  will  be  termed  the  differential  approxi- 
mation of  the  finite-difference  scheme,  while  an  expansion  up  to  second 
order  terms  in  time  and  space  is  termed  the  first  differential  approxima- 
tion [14,15,24]. 

The  stability  of  the  difference  schemes  may  be  investigated  by  means 
of  the  differential  approximation.  Such  investigations  were  made  by 
N.N.  Yaneriko  and  Y.I.  Shokin  for  one-dimensional  quasi-linear  equations  of 
the  hyperbolic  type  [24].  While  a strict  mathematical  foundation  has  not 
yet  been  supplied  for  the  case  of  non-linear  equations,  the  method  of 
differential  approximations  has  in  fact  been  used  here  [13]. 
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Taking  the  one -dimensional  case  for  simplicity,  let  us  describe  the 
first  differential  approximation  of  our  difference  scheme.  Take,  say 
ui+l’  write  as  faction  u(x  + Ax,t),  and  expand  each  term  of  the 
finite-difference  equations  in  Taylor  series  in  the  neighborhood  of  the 
point  (x,t) . 

For  instance,  in  computations  of  AM°  from  the  expressions  C4)  of 
second  order  accuracy  we  obtain 


dp  + 3pu 

at  ax 


0, 


apu  + a(P  + Pu_)  = . la  + ± (Pe  , 
at  ax  ax  ax  ax'  * 

or  when  using  egressions  (4)  of  first  order  accuracy 

M + ML  = 1 (£  lP.1, 
at  ax  ax  K ax'  * 

M . ± [uCp  . oE)  . - •£  *e  &&  * t|f)  C*#), 


(6) 


(7) 


ax 


ax 


ax 


where  e = |u|ax/2.  The  differential  approximations  may  be  written  down 
similarly  in  the  case  of  two-dimensional  problems. 

On  the  left-hand  sides  of  (6)  and  (7),  the  exact  expressions  of  the 
initial  differential  equations  have  been  obtained,  while  on  the  right  we 
have  the  terms  which  are  a consequence  of  "viscosity"  effects  in  the  dif- 
ference equations.  The  terms  involving  q result  from  the  explicit  intro- 
duction of  an  artificial  viscosity,  while  the  terms  involving  e are  due  to 
schematic  viscosity,  which  appears  when  the  exact  differential  equations 
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are  replaced  by  finite -difference  equations  ("expansion  errors") . 

It  may  easily  be  seen  that,  as  the  mesh  is  refined  (ax  -►  0) , e -*■  0 
and  the  equations  of  the  differential  approximation  reduce  to  the  exact 
set  of  initial  equations.  In  practical  computations  (due  to  Ax,  At,..., 
being  finite) , terms  containing  e always  appear  implicitly  in  the  differ- 
ence scheme  even  when  q = 0,  which  are  in  turn  analogous  to  the  dissipa- 
tive terms  of  the  Navier-Stokes  equations.  The  role  of  the  coefficient 
of  actual  viscosity  is  here  played  by  the  coefficient  e of  schematic  vis- 
cosity, which  depends  on  the  local  flow  velocity  and  the  size  of  the 
difference  mesh. 

In  the  two-dimensional  case,  it  follows  from  the  equations  of  momen- 
tum that  the  schematic  viscosity  (with  q = o)  has  the  tensor  form 

A = 7 

where  Ar  = Axi  + Ayj 
It  is  clear  from  (6')  that,  due  to  the  presence  of  the  vectors  Ar  and 
the  schematic  viscosity  does  not  possess  invariance  under  Galileo  trans- 
formations; in  practice  it  only  appears  in  zones  where  the  gradient  is 

large,  i.e.,  in  a shock  wave,  at  the  body  surface,  and  near  flow  separation. 
The  coefficient  of  schematic  viscosity  e (and  hence  the  width  of  the 
"smeared"  shock  wave)  then  depends  on  the  size  of  the  local  flow  velocity 
and  the  cell  size.  In  regions  of  smooth  flow,  where  the  gradients  of  the 
flow  parameters  are  relatively  small,  the  influence  of  the  schematic  vis- 
cosity is  negligible. 

It  will  be  shown  that  in  certain  cases  (when  expressions  (4)  of 


22 


the  first  order  accuracy  are  used  for  computing  AM0,  the  schematic  vis- 
cosity ensures  a stable  computation  even  in  the  absence  of  pseudo -viscosity 
q;  whereas  when  the  second-order  expressions  (4)  are  used  in  regions  where 
the  local  velocity  is  small  compared  with  the  velocity  of  sound,  the  intro- 
duction of  a term  with  q is  necessary  to  obtain  a stable  solution. 

2.4  Stability  of  the  Scheme. 

While  it  is  natural  for  different  types  of  difference  equations  to 
be  appropriate  at  various  stages , the  computations  become  strongly  unstable 
on  occasions,  and  rapidly  increasing  and  oscillating  solutions  appear, 
which  no  longer  reflect  the  behavior  of  the  solutions  of  the  initial  dif- 
ferential equations. 

The  difference  schemes  quoted  above  are  of  the  multi-layer  type, 
while  the  difference  equations  are  strongly  non-linear  with  variable  co- 
efficients. This  makes  it  impossible  to  employ  Fourier's  method,  devised 
for  linear  equations  with  constant  coefficients,  for  investigating  the 
stability  of  the  difference  scheme  as  a whole.  In  essence,  Fourier's  me- 
thod presupposes  that  the  equations  are  linearized  in  the  neighborhood  of 
the  flow  with  constant  parameters,  and  it  ignores  non-linear  effects  (in- 
fluence of  the  flow  gradients) , which  are  sometimes  the  true  sources  of 
the  instability. 

A heuristic  approach  will  therefore  be  employed  here  to  analyze  the 
stability  of  the  difference  schemes,  based  on  a consideration  of  their  dif- 
ferential approximations  [13,14,15]  and  appropriate  for  non-linear  equa- 
tions. In  this  approach,  we  determine  the  signs  of  the  coefficients 
("diffusion  coefficients")  in  the  dissipative  terms  of  the  differential 
approximation;  these  terms  contain  second  partial  derivatives  in  the  space 
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variables.  For  example,  a linear  equation  can  be  indicated,  such  that, 
when  the  value  of  the  coefficient  is  negative,  the  equation  of  the  differ- 
ential approximation  admits  a solution  which  is  exponentially  increasing 
in  time  (unstable)  [14]. 

In  short,  the  necessary  conditions  for  stability  are  obtained  here 
from  the  condition  o > 0 (parabolicity  condition) . In  the  case  of  linear 
equations,  the  results  of  a stability  analysis  obtained  by  means  of  the 
differential  approximation  are  exactly  the  same  as  that  obtained  by 
Fourier's  method. 

Let  us  examine  how  the  different  ways  of  writing  the  continuity  equa- 
tion (second  stage  of  the  computations)  contribute  to  the  instability, 
assuming  that  the  equations  of  momentum  and  energy  are  stable. 

If  AM0  is  determined  from  the  second  order  accurate  expressions  of 

(4) , we  find,  on  expanding  the  relevant  difference  equations  in  Taylor 

2 2 

series  and  retaining  terms  containing  3 p/3x  : 


3p  . 3pU  At  ( 2 2-.  3 p 

3t  + IF  - h - T (u  + c } • (8) 

If  AM°  is  evaluated  from  the  first  order  accurate  expressions  of  (4) , 
we  get 


3p  + 9pu  - 
3t  3x 


* 


+ c2) 


2 

AX  3u -| 

T~  3xJ 


f 


* 


(9) 


where  A^  and  A^  are  terms  of  the  first  differential  approximation  propor- 
tional to  ax  and  containing  the  first  derivatives.  In  our  case  [14,15], 
Ax  * 0.071;  At  * 0.0071;  pa  = 1;  u^  = 1 (10) 

In  practical  computations,  when  shock  waves,  contact  discontinuities  and 
rarefaction  wave  appear, 

pM  3 1;  |f£|  AX  < 0.3;  ||||  Ax  < 2. 
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It  follows  from  this  that  the  coefficient  of  3 p/3x  in  (9)  is  posi- 
tive, whereas  it  is  negative  in  (8),  i.e.,  scheme  (8)  is  computationally 
unstable,  while  scheme  (9)  is  stable. 

2.5  Some  Practical  Aspects  of  the  Method 

It  follows  from  the  very  character  of  the  construction  of  the  calcula- 
tion scheme  that  a complete  system  of  nonstationary  gas  dynamics  equations 
is  essentially  solved  here,  while  each  calculation  cycle  represents  a com- 
pleted process  in  calculating  a given  time  interval.  Besides,  all  the 
initial  nonstationary  equations  as  well  as  the  boundary  conditions  of  the 
problem  are  satisfied  and  the  real  fluid  flow  at  the  time  in  question  is 
determined. 

Thus,  the  "large  particles"  method  allows  us  to  obtain  the  properties 
of  nonstationary  gas  flows  and  as  a consequence  of  their  stability 
characteristics  their  asymptotic  state  as  well.  Such  an  approach  is 
especially  applicable  to  problems  in  which  a complete  or  partial  develop- 
ment of  physical  phenomena  with  respect  to  time  takes  place.  For  example , 
in  studying  transonic  gas  flows  and  flows  around  finite  bodies,  flow  in 
local  supersonic  zones,  and  separation  regions  develop  comparatively  slowly 
while  the  major  part  of  the  field  develops  rather  rapidly. 

In  contrast  to  the  FLIC  - method  [25]  our  investigation  is  wholly  de- 
voted to  systematic  calculations  of  a wide  class  of  compressible  flows 
involving  transonic  regimes;  discontinuity,  separation  and  "injected" 
flows,  etc. 

The  divergence  forms  of  the  differential  and  difference  equations  are 
considered  in  the  "large  particle"  method;  different  kinds  of  approximations 
are  used  in  the  1st  and  2nd  stages;  additional  density  calculations  are 
introduced  in  the  final  stage,  which  helps  to  remove  fluctuations  and  makes 
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it  possible  to  obtain  satisfactory  results  with  a relatively  small  network 
(usually  one  to  2 . 5 thousand  cells  are  used) . All  this  results  in  completely 
conservative  schemes,  i.e.,  laws  of  conservation  for  the  whole  mesh  region 
are  an  algebraic  consequence  of  the  difference  equations.  Fractional  cells 
are  introduced  for  the  calculation  of  bodies  with  a curvature  in  the  slope 
of  the  contour. 

The  investigation  of  these  schemes  (approximation  problems,  viscosity, 
stability,  etc.)  was  carried  out  for  the  zero,  the  first  and  the  second 
differential  approximations  [13-15].  These  investigations  show  that  the 
"large  particle"  method  yields  divergence-conservative  and  dissipative- 
-steady  schemes  for  "sweeping -through"  calculations. 

These  enable  us  to  carry  out  stable  calculations  for  a wide  class  of 
gas  dynamics  problems  without  introducing  explicit  terms  with  artificial 
viscosity.  It  may  be  of  particular  significance  in  studying  flows  around 
bodies  with  a curvature  in  the  slope  of  the  contour  since  the  ways  of 
introducing  explicit  terms  with  artificial  viscosity  are  different  for 
whole  and  fractional  cells.  Moreover,  by  varying  only  the  second  stage  of 
the  calculation  procedure  we  can  arrive  at  the  conservative  "particle- in-cell" 
method  so  that  the  calculational  algorithm  is  of  general  use. 

As  for  discontinuities  the  approximate  viscosity  in  the  schemes  (dissi- 
pative terms  in  difference  equations)  results  in  stable  calculations  with 
a "smearing"  of  shock  waves  over  several  computational  cells  and  the  forma- 
tion of  a wide  boundary  layer  near  the  body.  It  should  be  stressed  that 
the  magnitude  of  the  approximate  viscosity  is  proportional  to  a local  flow 
velocity  and  to  the  dimension  of  the  difference  mesh,  therefore  its  effect 
is  practically  evident  only  in  zones  with  large  gradients. 


26 


2.6  Results  of  Numerical  Computations 

Seme  computational  results  obtained  by  the  "large  particle"  method 
[20]  for  transonic  and  "supercritical"  flows  around  profiles,  plane  and 
axisynmetrical  bodies  will  now  be  described. 

For  purposes  of  this  discussion  the  supercritical  regimes  of  transonic  I 

flows  around  bodies  will  be  characterized  by  the  value  of  the  critical  Mach 
number  of  the  oncoming  flow  M^*  (i.e.,  when  a sonic  point  first  develops  on  / 

the  body)  as  well  as  by  the  extent  of  the  local  supersonic  zone  (as  com- 
pared to  a characteristic  dimension  of  the  body)  and  by  its  intensity  (say 
the  maximum  supersonic  velocity  realized  in  the  zone) . 

Figure  1 (series  l.a  - l.h)  presents  the  flow  field  patterns  (lines 
M = const.)  for  a 24%  circular  arc  profile  (v  = 0)  extending  from  purely 
subsonic  (Mm  = 0.6)  to  supersonic  regimes  (M^  = 1.5).  Successive  flow 

fields  for  increasing  depict  transition  through  the  critical  Mach  num- 

* 

ber  (here  M^  = 0.65),  and  the  formation  and  development  of  a local  super- 
sonic zone.  The  supercritical  flow  around  this  profile  is  observed  in 
Fig.  1 (b)-(g)  (0.65  < Mm  < 1).  One  can  distinctly  see  the  position  of  the 
shock  in  the  region  of  crowded  lines  M = const,  which  bound  the  local  super- 
sonic flow  together  with  the  sonic  line  (M  = 1) . The  region  of  subsonic 
velocities  is  located  behind  the  shock  wave.  When  the  velocity  of  the  on- 
coming flow  increases,  the  flow  disturbances  produced  by  the  body  die  out 
at  a large  distance  from  the  body.  With  Mm  > 0.9  the  zone  becomes  consi- 

a 

derable  both  in  size  and  in  intensity  (supersonic  velocities  up  to  M*  1.7 
to  1.8  are  attained)  and  in  case  of  sonic  flow  (Fig.  l.g)  the  sonic  lines 
(M  = 1)  extend  to  infinity. 

The  asymnetry  of  the  whole  flow  pattern  is  noticeable  (even  at  purely 
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subsonic  velocities  - Fig.  l.a)  which  results  from  non-potentiality  of  the 
flow  (super -critical  regimes)  and  from  the  presence  of  viscous  effects  as 
well  (formation  of  a wake  behind  the  body) . 

In  the  case  of  a supersonic  flow  around  a profile  (Fig.  l.h,  = 1.5) 
a shock  wave  ahead  of  the  body  develops  which  bounds  the  disturbed  region. 

Behind  the  wave,  subsonic  velocites  occur  in  the  vicinity  of  the  axis  of 
symmetry,  away  from  which  the  flow  velocity  along  the  contour  of  the  body 
increases  and,  as  a result,  a "terminal"  shock  occurs  near  the  stem  of 
the  body. 

Fbr  comparison  the  results  of  calculations  by  the  above  method  of  a 
flow  around  a 24%  axisymmetric  "spindle-like"  body  (v  = 1)  are  given  in 
Fig.  2 (2.a-2.h)  for  0.8  < < 2.5  . In  this  case  a critical  regime 

already  occurs  at  M*  = 0.86;  local  supersonic  zones  as  compared  to  the 
plane  case  are  less  developed  and  of  weaker  intensity  (for  example , values 
of  M ~ 1.3  to  1.4  are  realized),  although,  naturally,  the  main  features  of  a 
transonic  flow  are  quite  evident. 

In  Fig.  3 a comparison  is  given  between  the  flow  fields  calculated  by 
the  above  method  (solid  line)  and  those  of  the  Wood  and  Goodenm  experiment 
(dashed  line)  [26]  for  subcritical  (Fig.  3. a,  = 0.725)  and  supercritical 
(Fig.  3.b  = 0.761)  flows  around  a 121  profile  (results  of  both  the 

calculations  and  the  experiment  indicate  M£  = 0.74). 

The  analysis  of  internal  check  tests  as  well  as  the  results  of  com- 
parisons indicate  that  the  computational  error  of  the  "large  particle" 
method  usually  does  not  exceed  several  per  cent.  The  calculations  were 
carried  out  using  a Soviet  BESM-6  computer;  the  time  of  the  calculation 
in  this  case  did  not  exceed  an  hour. 

I 
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Figures  4 to  6 show  results  of  calculations  for  some  complicated  flows 
past  bodies  of  different  shapes  in  the  presence  of  discontinuities  in  the 
wake  as  well  as  under  the  influence  of  injection  of  a fluid  upstream  from 
the  front  surface  of  the  body.  Such  flows  are  of  great  practical  interest 
in  the  study  of  wakes  and  turbulence. 

The  results  of  numerical  experiments  for  the  investigation  of  flows 
with  injection,  are  given  in  Figs.  4a,  5 and  6b, c.  They  include  the  case 
of  the  interaction  of  a supersonic  flow  around  a finite  thick  circular  disk 
(Fig.  4a,  = 3.5),  a 241  body  of  revolution  (Fig.  5,  = 3.5),  and  a 

sphere  (Fig.  6b,  = 3.5;  Fig.  6c,  Mro  = 6)  with  a sonic  injection  stream 

(i.e.,  one  where  Mc  = 1.0;  pc  = 2.9,  uc  = 1.0,  vc  = 0)  issuing  upstream 
out  of  a nozzle  situated  on  the  axis  of  symmetry  of  the  body.  Fig.  6d 
presents  results  for  the  case  when  distributed  injection  of  the  flow  takes 
place  at  the  surface  of  a sphere.  In  all  the  figures,  streamlines,  shock 
waves,  horizontal  velocity  lines  (dots),  and  sonic  lines  (circles)  are 
indicated;  dashes  denote  lines  separating  the  main  flow  from  the  injected 
stream. 

The  action  of  the  jet  markedly  complicates  the  flow  pattern.  For 
instance,  in  the  flow  past  a cylinder  the  head  shock  wave  ABCD  (Fig.  4a)  is 
pushed  towards  the  oncoming  flow,  and  its  distance  from  the  body  increases 
significantly.  The  jet  issues  out  of  the  body  in  the  direction  of  axis  of 
symmetry  at  a sonic  velocity  and  expands,  forming  a local  supersonic  re- 
gion OLMNPO  which  is  closed  by  a triple -shock  intersection  (normal  front: 
ML;  oblique  front :MN;  transverse  front :MP) , having  a common  point  M.  A 
recirculation  zone  with  a complicated  vortex  structure  developes  in  front 
of  the  body;  sonic  line  BQ  is  situated  much  lower  in  comparison  to  flow 
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without  injection.  Behind  the  front  recirculation  zone,  a secondary  shock 
QC  is  formed  and  at  some  distance  from  the  body  merges  with  the  head  shock 
wave  ABCD  at  point  C. 

Behind  the  bodies  in  Figs.  4-6,  both  with  and  without  injection,  one 
can  observe  the  development  of  separated  zones  of  recirculating  flows.  In 
the  cases  considered,  these  zones  are  closed,  localized  in  the  wake  of  the 
body  and  separated  from  the  external  flow  by  a "non- flow"  line,  i.e. , a 
contact  surface  indicated  by  dashes  in  the  figures.  In  the  vicinity  of 
the  separation  (it  is  interesting  to  note  that  in  Fig.  4 the  separation 
point  is  situated  somewhat  lower  than  the  rear  shoulder  of  the  body)  a 
transverse  shock  wave  FF  develops.  Backward  recirculation  flows  are  essen- 
tially subsonic  and  rarefied  (gas  density  and  pressure  are  low) , so  that 
effects  of  viscosity  are  negligible. 

The  "large  particle"  method  has  also  been  applied  to  the  study  of  in- 
ternal gas  flows  and  diffraction  problems.  Fig.  7 presents  results  of 
computations  for  flow  through  a straight  channel  (v  = 0,  Fig.  7a)  and  a 
straight  tube  (v  = 1,  Fig.  7b)  in  the  presence  of  a central  body  (M^  =1.5) 
for  the  case  when  a triple  shock  intersection  is  formed  as  a result  of  the 
interaction  of  the  flow  with  the  upper  wall  (this  can  be  seen  by  an 
examination  of  the  lines  M = constant  in  Fi,*.  7a  and  rot  W = constant 
in  Fig.  7b). 

In  calculating  separated  flows,  the  cell  dimensions  of  the  "large 
particles"  were  changed  several  times  so  that  across  the  wake  of  a body  of 
size  R from  4 to  30  computational  intervals  were  used  (Fig.  8).  In  all 
cases  there  was  an  ample  reserve  of  computational  stability  (over  100  x 
Courant,  where  the  Courant  number  represents  the  ratio  of  the  time  step  to 
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the  space  width  of  the  cell) . 

Fig.  8 shows  base  flows  behind  an  axisymmetric  cylinder  (M^  = 2.0)  for 
R = 14 Ay.  A gradual  development  of  the  flow  in  time  is  shown  (in  dimen- 
sionless units)  from  tn  = 21  to  tn  = 31,  when  the  zone  is  practically  loca- 
ted. Streamlines  are  represented  by  solid  lines;  velocity  vectors  by 
arrows.  It  follows  from  this  diagram  that  at  tn  = 25  the  flow  has  already 
been  formed  but  still  continues  to  "breathe".  It  is  interesting  to  note 
that  similar  flow  patterns  were  obtained  with  denser  meshes  (which  is  quite 
important)  and  the  zone  ,fbreathing"  i.e.,  changes  of  its  dimensions,  inter- 
nal structure,  and  other  features  of  the  flow  occurred  approximately  at 
the  same  time  intervals,  tn,  in  the  various  approximations. 

The  development  of  flow  separation  in  the  case  of  strong  interaction 
seems  to  be  explained  by  the  fact  that,  as  a result  of  viscosity  (compu- 
tational) effects  and  the  treatment  of  the  boundary  conditions,  close  to 
no-slip  conditions  are  realized  on  the  body  itself;  a fairly  thick  boundary 
layer  forms  around  the  body  surface  (comparable  to  the  width  of  the 
body  at  its  tail) , and  this  layer  then  separates  from  the  body  surface 
and  forms  a near  wake  flow  with  complicated  vertical  structure  behind  the 
base  of  the  body.  It  must  be  empahsized  that,  while  the  boundary  layer  is 
in  fact  the  result  of  viscosity  effects  in  the  scheme,  in  the  wake  itself 
the  influence  of  the  approximation  viscosity  e ( which  is  proportional  to 
the  local  velocity  and  the  size  of  the  computational  mesh,  see  above)  is 
quite  small,  since  in  these  zones  only  small  values  of  the  subsonic  velo- 
cities are  realized,  while  computations  with  different  size  meshes  revealed 
only  a slight  change  (within  the  limits  of  one  step)  in  the  zone  contour. 

The  fact  that  the  solution  does  not  strongly  depend  on  viscosity 
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e - p/u/h)  shows,  by  the  way,  that  flows  corresponding  to  high  Reynolds 
numbers  can  be  treated  by  our  methods  of  analysis.  Thus,  our  calculations 


of  separated  zones  might  give  quantitative  information  for  "limiting" 
flows  (Re  -*■  »)  as,  for  example,  the  calculation  of  shock  waves  by  a scheme 
including  viscous  effects.  Naturally,  the  accuracy  of  determining  the 
characteristic  features  of  such  zones  can  be  further  increased,  if  necessary 
by  using  the  results  of  preliminary  calculations  (e.g.,  the  position  of 
separation  and  closure  points,  a zone  contour,  etc.)  as  initial  data. 

However,  it  should  be  pointed  out  that  in  the  calculations,  the 
flow  parameters  on  the  front  part  of  the  body  are  determined  com- 
paratively quickly,  while  local  supersonic  zones  and  separation  regions 
continue,  as  mentioned  above,  to  "breathe".  This  may  be  due  to  the  physi- 
cal (non- stationary)  character  of  the  phenomenon  itself.  Hie  application 
of  the  difference  scheme  of  calculations  prescribed  by  the  non-stationary 
method  of  "large  particles"  appears  to  be  especially  well  suited  for  such 
a case. 
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3.  Computation  of  Incompressible  Viscous  Flows. 

At  the  present  time,  a fairly  large  number  of  numerical  methods  of 
solving  the  Navier-Stokes  equations  (governing  viscous  incompressible  flows) 
are  known.  Most  of  these  methods  were  developed  for  equations  expressed  in 
terms  of  the  stream  function  y and  vorticity  a>. 

A common  disadvantage  of  these  methods  is  the  need  for  some  form  of 
a boundary  condition  (the  Tom  condition)  for  the  vorticity  on  the  solid  surface, 
which  is  absent  in  the  physical  formulation  of  the  problem.  The  rate  of 
convergence  of  numerical  algorithms  is  limited  by  the  presence  of  an  addi- 
tional iteration  imposed  by  this  boundary  condition  on  the  surface  vor- 
ticity. 

Moreover,  the  obvious  limitation  of  methods  of  solving  the  (y,ai) -system, 
connected  with  their  inapplicability  for  cases  of  three-dimensional  viscous 
flows  and  compressible  gas  flows,  accounts  for  the  recent  interest  in  the 
numerical  solution  of  Navier-Stokes  equations  expressed  in  natural  variables: 

||  + $-v)$  = -flp  + \>A$,  = 0 (11) 

where  p - pressure,  ^ - vector  of  velocity,  v - coefficient  of  kinematic 
viscosity. 

Using  the  main  principles  of  the  "large  particle"  method  V.A.  Gushchin 
and  V.V.  Shchennikov  [16,17,29]  studied  viscous  incompressible  gas  flows 
with  variables  "velocity-pressure"  by  means  of  a numerical  scheme  of 
splitting  analogous  to  the  SMAC  method  [27]. 

3.1  Description  of  the  Splitting  Method 

Let  us  now  consider  the  scheme  of  difference  approximations  of  equa- 
tions (11)  which  enable  calculations  by  a single  algorithm  for  plane, 
axisymmetric  and  three-dimensional  flows  of  a viscous  incompressible  fluid. 
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For  this  purpose  consider  the  following  scheme  of  splitting  a time  cycle: 
Stage  I - determination  of  intermediate  values  of  velocities 

(^  - ^n)  /t  = -C^n  • v)^n  + vA^n  ; (12.1) 

u Stage  II  - calculation  of  the  pressure  field 

sp  . 4»/t  , (v  • ?n*1 . rf1*1  . 0)  ; (12-2) 

Stage  III  - determination  of  the  final  values  of  velocities 

$n+1  = $ - x*ty.  (12.3) 

Stages  I and  III  lead  to  the  realization  of  the  Navier-Stokes  equation 
and  stages  II  and  III  are  the  conditions  of  solenoidality  (second  equation  of 
(11)).  Consequently  at  stage  I the  evolution  of  the  velocity  field  is  ac- 
complished only  by  convection  and  diffusion  so  that  the  resulting  ^-field 
does  not  satisfy  the  continuity  equation  (i.e.,  D + 0).  Therefore  it  is 
necessary  to  change  ("to  correct")  the  field  ^ at  the  expense  of  a pres- 
sure gradient  p so  that  = 0 (stage  III)  where  p is  found  by  solving 
the  Poisson  equation  (stage  II). 

For  a proportional  calculational  mesh  a two-dimensional  difference 
scheme  of  second  order  accuracy  with  respect  to  space  is  presented  in 
[17].  The  main  difficulties  of  the  numerical  realization  of  the  scheme 
involve  the  calculation  of  a pressure  field  and  the  formulation  of  boundary 
conditions. 

It  should  be  noted  that  in  some  works  boundary  conditions  at  a solid 
surface  are  replaced  by  the  projection  of  an  equation  of  motion  onto  the 
normal  to  the  surface  at  the  boundary  points.  This  substitution  reduces  the 
efficiency  of  the  numerical  methods  since  these  conditions  are  not  available 
in  the  physical  formulation  of  the  problem. 

In  [28]  there  is  proposed  an  original  modification  of  boundary  conditions 
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in  the  MAC  method,  which  enables  homogeneous  boundary  conditions  to  be 
provided  for  pressure.  Moreover,  in  the  SMAC  method  [27]  and  the  modified 
MAC  method  [28],  due  to  the  difference  schemes  chosen,  the  realization  of  the 
condition  of  no-slip  necessarily  results  in  the  determination  of  the  vor- 
ticity  value  on  a solid  surface  which  satisfies  the  Tom  condition  to  first 
order  accuracy.  In  addition,  the  no-slip  condition  in  the  SMAC  method  does 
not  provide  a balance  of  forces  on  a solid  surface.  The  error  in  this 
case  is  of  the  order  of  0(v). 

An  essential  point  of  the  method  proposed  is  the  choice  of  boundary 
conditions.  From  the  viewpoint  of  the  solution  of  problems  of  a viscous 
incompressible  flow  around  bodies  of  finite  dimensions  we  can  distinguish 
two  basic  types  of  boundary  conditions:  conditions  on  a solid  surface  and 
those  on  a line  sufficiently  remote  from  the  body.  Let  us  consider  each 
of  these  conditions  in  further  detail. 


Boundary  conditions  on  a solid  surface: 


"U  ■ 0 

'i+.-h ' 0 

from  the  latter  it  follows 


(nonpenetration  condition) , 
(no -slip  condition) ; 


Vh.o  ■ % * ^ ♦ °0i3) 


Condition  (14)  makes  it  possible  to  determine  the  boundary  value  for  u 
of  second  order  accuracy  with  respect  to  internal  field  points.  This  avoids 


the  alternative  of  introducing  a layer  of  fictitious  cells  (inside  a solid 
body),  which  in  schemes  of  the  MAC,  SMAC  and  modified  MAC  types  [23]  gives 


rise  to  only  first  order  accurate  values  of  surface  vorticity.  Note  that  in 
the  limits  of  the  approach  proposed  it  is  unnecessary  to  calculate  the 
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vorticity  on  the  solid  surface.  The  latter  can  be  determined  from  a calcu- 
lated velocity  field  using  some  of  the  difference  representations  of  the 

vorticity 

0)  = 3u/3y  - 3v/3x 
at  boundary  points. 

Boundary  conditions  on  a line  remote  from  the  body  represent  an  undis- 
turbed flow;  for  UJI  OX  this  has  the  form 

vi,N+*5  = °»  ui+3s,N  = U~ 

In  calculating  the  pressure  field,  homogeneous  boundary  conditions 
are  attained  following  the  approach  of  [28].  Corresponding  to  v?+^,  = 0 
(for  the  case  of  a solid  surface)  and  ^ = 0 (for  the  case  of  a line 
remote  from  a body)  we  have  from  the  finite -difference  approximation  of 
(12.3) 

vi,-h  = E (pi,0  ‘ Pi,-^,  vi,N+*s  = E (pi,N*l  ' Pi.N5*  fl5) 

Taking  account  of  (15)  it  is  not  difficult  to  write  down  now  a difference 
equation  for  calculating  the  pressure  at  boundary  cells  [17]. 

The  stationary  solution  of  the  system  of  equations  (12)  is  obtained  as 
a result  of  repetition  of  the  above  stages  until  the  following  criterion 
is  fulfilled 

.n+1 


±9J 


The  stability  can  be  investigated  stage-by-stage.  A stability  cri- 
terion for  the  first  stage  can  be  obtained  from  the  first  differential 
approximation  (condition  of  a-parabolicity) . With  regard  to  equations 
(12a)  the  first  differential  approximation  is  [17]: 
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The  stability  criterion  of  the  difference  scheme  employed  follows  from  (16) 
t < 4v/(u ^ + v2) 

Eliminating  p from  (12.2)  and  (12.3),  the  unconditional  stability  of  the 
second  and  third  stages  is  easily  demonstrated  by  Fourier's  method. 

Thus,  the  proposed  difference  scheme  enables  us  to  calculate  a flow 
without  prescribing  vorticity  and  pressure  at  the  body  surface.  This 
markedly  increases  the  accuracy  of  the  calculations.  Results  of  the  cal- 
culations attest  to  its  effectiveness.  This  difference  scheme  (of  second- 
-order  accuracy)  provides  a single  algorithm  for  calculating  viscous  in- 
compressible flows  around  plane,  axy symmetrical  and  three-dimensional  bodies 
of  complex  configuration  as  well  as  internal  flows  in  a wide  range  of 
Reynolds  nunbers  [17]. 

3.2  Results  of  Numerical  Computations 

Solutions  of  a whole  series  of  problems  of  external  hydrodynamics  were 
obtained  by  the  method  described.  Viscous  incompressible  flows  over  a 
wide  range  of  Reynolds  nunbers  (1  < Re  < 10^)  were  studied  around  different 
bodies  of  finite  dimensions:  a rectangular  slab  and  a cylinder  of  finite 
length  whose  axis  is  parallel  to  the  free  stream  [29],  a 
sphere  and  a cylinder  with  the  axis  perpendicular  to  TT^,  a rectangular 
parallelepiped  (a  three-dimensional  flow)  [7],  as  well  as  bodies  of  more 
complex  form. 

Fig.  9 shows  the  steady-state  streamline  patterns  around  a circular 
cylinder  (two-dimensional  flow)  for  Re  “ 1,10,30  and  50 (Re  - 2Rvw/v,  vfoere 
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R is  the  cylinder  radius).  The  streamline  patterns  in  flow  around  a cylin- 
der for  Re  = 103  are  shown  in  Fig.  10  at  the  instants  tj  = 162,  t2  = 166  and 
tj  = 170,  respectively.  In  the  last  case  a non-steady  flow  pattern  is  ob- 
served (there  is  a definite  growth  of  the  stagnation  zone  and  at  some  in- 
stant of  time  occurs  a "collapse"  and  ejection  of  fluid  from  the  stagnation 
zone).  This  result  probably  requires  further  study. 

Figs.  11  and  12  present  results  for  unsteady  three-dimensional 
incompressible  viscous  flow  around  a cube  (of  dimension  2a)  when 

the  oncoming  flow  velocity  is  parallel  to  an  axis  OX.  Due  to  the  pre- 
sence of  two  planes  of  symmetry  (OXY  and  OZX)  the  calculation  is  carried 
out  only  in  the  positive  quadrant  OXYZ  (Fig.  11) . The  properties  of  the 

flow  are  illustrated  by  the  velocity  profiles  u (parallel  to  the  vector  U ) 

00 

at  various  cross-sections  Q (x  = const.).  Fig.  11  shows  the  velocity 
profile  u in  the  undistrubed  flow  (x  = -«)  and  for  a section  x = 3a  with 
Re  = 1 (Re  = 2a  / v) . Fig.  12a  illustrates  for  various  Re  the  spatial 
change  of  the  velocity  profiles  u at  several  sections  downstream  of  the 
body  (section  coincides  with  the  rear  face  of  the  cube  x = 2a;  the 
distance  between  the  sections  is  constant.  Ax = 0.5a).  Fig.  12b  shows  the 
evolution  with  time  (1.0  < t < 1.29)  of  the  velocity  field  for  the  section 
x = 4a.  It  follows  from  Fig.  12,  in  particular,  that  with  Re  = 40  and  100 
a reverse -circulation  zone  (u<  0)  develops,  and  after  an  elapse  of  time  a cer- 
tain flow  stabilization  is  observed.  The  reader  is  referred  to  [16,17,29] 
for  further  details  concerning  results  of  this  calculation. 
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L Flow  of  a Viscous  Compressible  Gas  (Conservative  "Flow  Method") 

The  calculation  of  viscous  conpressible  gas  flows  was  performed  by 
L.I.  Severinov  and  A. I.  Babakov  with  the  approximation  of  conservation 
laws  represented  in  integral  form  for  each  cell  of  the  calculation  scheme 
("flow"  method)  [18].  Conservation  laws  for  mass,  momentum  and  energy  of 
a finite  volume  have  the  form: 

^=-//$FaS,  F = {M,X,Y,Z,E}  , (17 ) 

sn 

where  sfi  - is  the  lateral  surface  of  volume  cell  ft;  M - mass,  X,Y,Z,  - 
momentum  components  and  E - energy  terms  in  ft,  respectively,  and  $p  is  a 
flow  density  vector  for  each  of  the  quantities.  Eqs.  (17)  take  account  of 
the  boundary  conditions  and  are  solved  numerically  for  each  cell  of  the 
computational  domain. 

If  the  values  of  M11  = M(tn),  are  known  at  the  instant 

tn  = nt,  where  t is  a time  integration  step,  at  time  tn+1  = (n+1)  t 

2 

these  quantities  can  be  calculated  with  error  0(t  ) as  follows  [18]: 

Fn+1  = F°  - t //  • Zs  . (18) 

Sft 

Supplementary  conditions,  the  form  of  which  depends  on  the  problems 
posed,  must  make  it  possible  to  determine  the  flow-density  vectors  on  the 
boundary  of  the  domain  in  which  the  solution  is  sought.  The  three-dimen- 
sional coordinate  system,  the  shape  of  the  cell  volumes  ft,  the  methods  of 
determining  the  field  variables  and  their  first  derivatives  on  the  surfaces 
sQ  must  be  chosen  in  such  a way  as  to  ensure  stability  and  monotonicity  of 
the  difference  method,  as  well  as  a fairly  simple  approximation  of  the  in- 
tegrals in  the  system  (18). 
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In  the  solution  of  a specific  problem  the  integrals  in  (18)  are  cal- 
culated on  separate  segments  of  the  surface  which  are  the  boundaries  be- 
tween two  adjacent  volumes  n.  Depending  on  the  directions  of  the  flow 
density  vectors  the  values  of  M,X,Y,Z  and  E vary  (they  increase  in  some 
cells  and  decrease  in  others)  by  quantities  determined  by  the  flows  of  mass, 
momentum  and  total  energy  through  the  corresponding  segments  of  the  boundary. 
Apart  from  round-off  errors  this  calculation  method  cannot  lead  to  the  loss 
or  generation  of  the  quantities  (M,X,Y,Z,E)  due  to  computational  errors. 
Therefore,  the  flow  method  is  conservative  with  respect  to  mass,  momentum 
components  and  total  energy  [18]. 

In  the  finite-difference  approximation  of  Equation  (18)  Stokes'  assump- 
tion of  the  equality  of  the  mean  values  of  the  principal  stresses  (with 
a reversed  sign)  and  pressure  has  been  used.  If  the  field  variables  are 
sufficiently  smooth  and  the  assumptions  used  in  calculating  the  mass,  mo- 
mentum and  energy  flow  density  vectors  are  satisfied,  the  conservation  laws 
(17)  imply  the  complete  Navier-Stokes  equations  for  a compressible  gas, 
if  the  volume  ft  is  arbitrary. 

4.1  Computational  Problems 

We  will  now  consider  some  problems  involved  with  the  numerical  inves- 
tigation of  equation  (17) . 

Knowing  the  values  of  M,X,Y,Z,E  for  each  cell,  we  can  calculate  for  a 
given  small  volume  cell  fi  fixed  in  space  the  average  densities  (in  each 
cell)  of  the  given  quantities  p,£,n,c,e  : 

p = M/n,  £ = x/n,  n = Y/q,  x,  = z/n,  e = E/n  . 

From  these  functions  it  is  easy  to  arrive  at  generally  accepted  field 
varibles,  i.e.,  the  components  u,v,w  of  the  velocity  f and  the  specific 
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internal  energy  of  the  gas  e: 

u = C/p,  v = n/p,  w « e/p,  e = e/p  - (u2  + v2  + w2)/2. 

Applying  certain  procedures  of  interpolation  and  numerical  differentiation,  we 

determine  the  value  of  the  field  variables  and  the  first  derivatives  of 


u,v,w,e  on  the  boundaries  s of  the  cells  n [18]. 

In  determining  the  values  of  all  the  function  [except  the  distribution 
densities  p,C,n,e,e)  and  their  first  derivatives,  symnetric  formulas  were 
used,  for  example. 


W (Vi,k  * ViP  / 2- 


C— ) 

axm^,k 


3u 

^m+Js  k = (l<m*k+1  ’ Vk-1  + Um+1  ,k+l  ' Vl.k-P  1 4h2  * 

However  asymmetric  formulas  were  used  to  calculate  the  density  values  of  the 
distributions  p ,£ ,n , e , e : 

. ‘ '■Vl,k  if  um+^,n  » 0 < 

‘Wi,k  'l-V,k'“-Vk  if  ‘ ° ' 

These  approximations  ensure  second-order  accuracy. 

In  approximating  flow  density  vectors  ($p  an  essential  aspect  of  the 
method  required  that  the  distribution  densities  of  additive  characteristics 
such  as  densities  F are  calculated  on  the  boundary  sfi  of  volume  fi  in  a 
non- symmetrical  way  (extrapolation  in  the  direction  of  gas  flow) ; while  the 
other  parameters,  e.g.,  pressure  and  transfer  velocities  u,v  of  additive 
characteristics  are  calculated  according  to  symmetric  formulas  in  the  vis- 
cous stress  tensor  and  in  the  thermal  conduction  law.  We  believe  that  this 
treatment  accounts  for  the  region  of  influence,  which  is  an  important  factor 
in  the  investigation  of  complex  physical  flow  patterns.  The  presence  of  a 
"convective"  transfer  gives  unequal  weighting  to  the  space  directions  and  it 
is  desirable  to  take  account  of  this  fact  in  constructing  the  difference 
scheme. 

The  integral  form  of  the  conservation  laws  essentially  requires  the 
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approximation  of  derivatives  of  one  order  lower  than  that  required  for 
numerical  solutions  of  the  Navier-Stokes  equations.  By  its  formulation  the 
"flow"  method  is  conservative  with  respect  to  mass,  momentum  and  total 
energy,  both  locally  (for  each  cell  of  a difference  mesh!)  and  integrally, 
i.e.,  for  the  whole  computational  space  [18].  As  follows  from  (17),  con- 
servativeness results  from  the  fact  that  this  approach  is  based  upon  the 
difference  approximation  of  conservation  laws  written  down  for  each  cell  in 
terms  of  surface  integrals  of  vectors  of  flow  densities  ($p,  i.e.,  the  con- 
servation law  which  governs  an  arbitrary  gas  volume.  Indeed,  while  solving 
a given  problem,  surface  integrals  in  (17)  are  calculated  on  separate  sur- 
face segments  s^  which  constitute  boundaries  between  two  adjacent  volumes  Q. 
Depending  upon  the  direction  of  flow  vectors  the  values  F = (M,X,Y,E)  vary 
(they  increase  in  some  cells  and  decrease  in  others)  and  acquire  new  values 
determined  by  flows  of  mass,  momentum  and  total  energy  across  coninon  bounda- 
ries of  the  cells.  . 

It  should  be  noted  that  the  "flow"  method  is  essentially  another  de- 
velopmei  the  "large  particle"  method.  Ihe  difference  formulas  of  the 
"flow"  meti  .od  can  be  deduced  by  using  the  splitting  scheme  (3)  - (5)  for 
the  "transfer"  of  the  components  of  quantities  F. 

4.2  Examples  of  Calculations 

The  "flow"  method  has  been  applied  to  the  systematic  study  of  the 
characteristics  of  viscous  compressible  gas  flow  around  bodies  of  finite 
dimensions  over  a wide  range  of  Reynolds  numbers  Re.  Although  the  method 
formally  "works"  even  with  large  values  of  Re,  the  results  are  reliable 
only  when  the  boundary  layer  thickness  is  much  greater  than  the  size  of  the 
computational  mesh.  It  should  be  emphasized  that  the  division  of  the 
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vector  of  flow  density  ($p  into  convective  and  viscous  components  facilitates 
application  of  the  algorithm  to  the  calculation  of  an  ideal  gas  flow,  as 
well. 

The  results  given  below  represent  the  asymptotic  steady  state  obtained 
by  the  method  for  the  solution  of  a stationary  boundary  value  problem.  The 
investigation  of  a linear  model  reinforced  by  the  results  of  the  calculations 
themselves  showed  that  the  second-order  accurate  difference  scheme  exhibited 
conventionally  stability  and  monotonicity  [18].  The  reliability  of  the 
numerical  results  is  examined  for  a general  case  by  subdividing  the  compu- 
tational cell,  by  using  various  forms  of  boundary  conditions,  by  comparing 
with  the  results  of  other  calculations  as  well  as  with  measurements  taken 
from  an  experiment  [18,30]. 

Some  details  of  the  flow  past  a sphere  (separation  zones  of  a reverse-cir- 

4 

culation)  at  Mro=20  and  550  < Re^  < 10  are  given  in  Figs.  13  and  14.  Fig.  14 
shows  the  behavior  of  lines  p = const,  in  a separation  zone  behind  the  sphere 
with  Re  = 10^  and  1500  (Re  = Rv  /v). 

OO  OO  00 

Fig.  15  illustrates  the  variation  of  density  across  the  shock  layer  in 
the  neighborhood  of  the  forward  stagnation  point  (x  = 3°)  for  75  < Rero  < 10^ 
and  compared  with  that  of  an  ideal  gas  (Mm  = 20,  k = 1.4,  Rem  = °°) . With 
increasing  Re  the  density  in  the  shock  layer  approaches  its  limiting  value 
in  a viscous  thermally  non-conducting  gas;  the  tendency  towards  the  forma- 
tion of  a shock  wave  is  distinctly  seen. 

The  pressure  distribution  along  a blunt  body  (relative  to  the  pressure 
at  the  stagnation  point,  x = 0)  is  shown  in  Fig.  16:  "lines"  designate  the 
results  obtained  by  the  "flow"  method  (M^  = 6.05,  Re^  = 6.43*10^),  "crosses" 

- experimental  data  (G.M.  Riabinkova)  and  "circles"  - the  results  for  an 
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ideal  gas  (O.M.  Belotserkovskii  [31]).  Very  good  agreement  is  observed 
between  the  data.  In  this  way  transition  from  the  "viscous"  equations  (17) 
to  the  limit  of  an  ideal  gas  is  obtained. 

The  results  of  calculations  show  that  the  "flow"  method  makes  it  pos- 
sible to  study  viscous  compressible  gas  flows  around  finite  bodies  over  a 
wide  range  of  flow  regimes  (including  separation  zones)  up  to  large  Rey- 
nolds numbers. 


. Statistical  Model  for  the  Investigation  of  Rarefied  Gas  Flows 

A statistical  variant  of  the  method  of  large  particles  has  been 
investigated  by  V.E.  Yanitsky  [21-23]  for  the  solution  of  the  Boltzmann 
equation.  The  main  difficulty  of  this  problem  is  the  development  of  a model  of 
behavior  of  the  gas  medium  consisting  of  a finite  number  of  particles.  This 
study  combines  the  splitting  of  the  "large  particle"  method  in  terms  of 
Bird's  statistical  treatment  [32,33]  with  Kats'  ideas  [34]  about  the  exis- 
tence of  asymptotically  equivalent  models  to  the  Boltzmann  equation. 

As  is  typical  of  "particle- in-cell"  methods,  the  medium  is  simulated 
by  a system  containing  a finite  number  N of  particles  of  fixed  mass.  At  a 
given  instant  of  time  t each  cell  j contains  N(a,j)  particles  endowed  with 
certain  velocities.  The  main  calculation  cycle  comprises  two  stages: 

- at  the  first  stage  particles  only  collide  with  their  counterparts 
in  a cell  (collisional  relaxation]  and 

- at  the  second  stage  they  are  only  displaced  and  interact  with  the 
boundary  of  a reference  volume  and  with  the  surface  of  a body(collisionless 
relaxation) . 

The  main  distinction  between  the  model  suggested  in  [21-23]  and  Bird's 
model  lies  in  the  fact  that  at  the  first  stage  of  the  calculation  each  group  of 
N particles  in  a cell  is  regarded  as  Kats'  statistical  model  for  an  ideal 
monoatomic  gas  consisting  of  a finite  number  of  particles  in  a homogeneous 
coordinate  space.  Molecular  collisions  are  calculated  by  the  Monte-Carlo 
method  from  the  main  equation  of  Kats'  model,  which  correctly  determines  the 
time  between  particle  collisions  in  accordance  with  collision  statistics 
for  an  ideal  gas. 

In  contrast  to  previously  reported  versions  of  Bird's  method  [32,33] 
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the  approach  in  question  [21-23]  is  a rigorously  Markovian  process.  The 
main  equation  of  this  formulation  is  linear  (unlike  the  Boltzmann  equation) , 
which  substantially  simplifies  numerical  realization  of  the  algorithm.  The 
feature  of  molecular  chaotic  motion  implies  that  Rats'  model  is  asymptoti- 
cally equivalent  to  the  Boltzmann  equation  without  the  convective  derivative. 
Integration  of  Rats'  main  equation  results  (with  accuracy  up  to  the  reali- 
zation of  the  assumption  of  molecular  chaotic  motion)  in  the  Boltzmann 
equation. 

In  the  realization  of  the  second  stage  of  the  calculation  for  the  dis- 
placement of  the  particles,  the  numerical  algorithms  [21-23]  employ  incomplete 
information  about  the  position  of  particles  in  coordinate  space.  This 
reduces  storage  requirements  in  the  processor  memory,  which  significantly 
increases  the  speed  of  the  computations.  The  method  can  just  as  well  be 
realized  in  a two-  or  three-dimensional  coordinate  space. 

Let  us  review  here  the  principal  aspects  of  the  statistical  "particle- 
-in-cell  method  [21-23]. 

5.1  Description  of  the  Statistical  "Particle- in-Cell"  Method. 

We  suppose  that  the  problem  of  a rarefied  gas  flow  around  a body  can 
be  solved  by  means  of  a distribution  function  with  a monoatomic  gas.  Then 
any  macroparameter  of  a gas  flow  t(t,x)  related  to  a molecular  feature 
M'(c)  is  a functional  of  the  form 

v(t,x)  = /¥  (c)-f(t,x,c)  dc  , 

n(t,x) 

where  f(t,x,c)  is  a molecular  distribution  function  in  a 6-dimensional  space 
(x,c)  of  the  coordinates  and  velocites  of  the  particles. 

If  n denotes  the  region  of  a control  volume  and  r - the  boundary  of  n 
which  emcompasses  a body  surface  as  well,  then  the  problem  is  reduced  to 
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that  of  obtaining  the  solution  of  Boltzmann's  equation 


||  + c si  = /(f  . f!  - 
ax  3 

ffj)  gda  • dCj  , 

(19) 

subject  to  the  initial  parameters 

f(t  + 0,x,c)  = fQ(x,c), 

xen , -<»  < c < + oo 

x,y,z 

(20) 

and  boundary  conditions 

f(t,xr,c)  = /K(c,c1)f(t,x?,c] 

Hep  cn(xr)>0,  CjSfXj.)  < 0 . 

(20') 

Here  nC^,)  is  the  normal  to  the  surface  r at  point  xLer  directed  toward 
the  interior  of  volume  fi;  the  kemal  shape  k is  derived  from  the  interaction 


law  of  the  "gas -surface". 

In  deducing  Boltzmann's  equation  the  following  assumptions  are  made. 

1)  mechanics  of  collisions  are  described  in  a classical  way; 

2)  force  fields  of  molecules  are  spherically  symnetric; 

3)  only  binary  collisions  are  considered  (two  molecules  take  part  in 
any  collision) ; 

4)  molecules  move  randomly  (the  hypothesis  of  molecular  chaos  is  valid, 
i.e.,  the  distribution  function  of  molecular  pairs  f2(t,x,CpC2)  = 

= f^(t,x,Cj)  • fj(t,x,c2)  which  implies  statistical  independence 
of  particles); 

5)  the  collision  time  is  negligibly  small. 

The  difficulty  in  constructing  the  solution  of  the  Boltzmann  equation 
in  nonlinear  integro-differential  form  results  both  from  the  large  nunber 
of  independent  variables  (there  are  seven  of  than  in  the  general  case:  time, 
geometrical  coordinates  and  molecular  velocity  components)  and  from  the  com- 
plex structure  of  the  integral  of  the  collisions.  Quadratic  nonlinearity  in 
the  integrands,  their  strong  dependence  upon  distribution  functions  (deter- 
mined by  the  values  of  molecular  velocities  after  a collision)  as  well  as 
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the  high  level  of  multiplicity  of  integration  (equal  to  five  in  a general 
case)  and  the  complex  formulation  of  the  boundary  condition  (20)  - these  are 
the  main  features  which  complicate  a direct  solution  of  the  Boltzmann 
equation  (19)  and  the  application  of  ordinary  nunerical  algorithms. 

For  an  approximate  solution  of  the  problems  formulated  in  this  fashion 
we  shall  construct  a statistical  model  of  an  ideal  monotomic  gas  consisting 
of  N particles*  with  coordinates  r^,  and  velocities  c^(i  = 1,2,...,N)  so 
that  the  equation  of  evolution  of  the  model  approximates  equation  (19) , the 
only  additional  assumption  being  that  of  molecular  chaos: 

f2(t,x,c1,c2)  = ^(t.x.Cjh  f-^t.x,^)  , (21) 


where 


fc(t,x,c1,...,cs)  - (N-g)l  Fs(t,r1,...,rs,c1,...,cJ 


s'  ”1’  ’S'  IN-SJ!  -S'-’-l’  - ’ S'  1’  * S' 

and  with  ^ = r2  = ...  = rg  = x,  Fg  being  a s-partial  function  of  distribu- 
tion in  a phase  space  of  6N  dimensions. 

If  {$(t),c(t)>  = {^(t)  ,cx  (t)  ; . . . ; rN(t),cN(t)}  designates  the  model 
state  at  time  t,  the  problem  solution  is  then  reduced  to  a numerical  reali- 
zation of  a finite  number  of  trajectories  {ft(t),c(t)}  with  initial  parame- 
ters corresponding  to  (20) ; the  modelling  of  particle  interaction  with  the 
boundary  r is  accomplished  in  accordance  with  the  given  kernel  * of  (20’). 
Having  calculated  a nunber  of  trajectories  one  can  obtain  any  macroparameter 
using  adequate  estimates  of  the  integrals  by  the  Monte  Carlo  method. 


* A real  gas  is  modelled  by  an  ensemble  of  about  a thousand  rigid  sphere-like 

molecules  that  can  be  regarded  as  typical  representatives  of  many  trillions 
12 

(10  ) of  molecules,  e.g.,  in  the  study  of  phenomena  occur ing  in  a real 
shock  wave  [22]. 
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The  synthesis  of  the  basic  ideas  of  splitting,  the  "particle"  method 
and  Rats'  statistical  model  enable  us  to  construct  the  desired  model 
($(t),c(t)}  for  a space  - inhomogeneous  case  when  3f/3x  t 0. 

Let  us  suppose  that  at  time  interval  t^fa  = 0,1...)  in  a cell  with 
center  Xj  ( j = l,2,...,j)  there  are  N(a,j)  particles  with  velocities 
<?l»--*»*n(a  j) ^ center  xj  a ceH  in  which  a particular  particle 
is  situated  is  taken  as  a coordinate  r^  of  a particle  i.  The  state  of  such 
a gas  model  {ft,c}  is  uniquely  defined  by  a sequence  of  J points  of  the  form: 
($(t),c(t)}  ~ {N(a,j);  Cs,...,CNkj)>  , j = 1,2,. ..,J  , 

J 

N = E N(a,j) 

The  evolution  of  this  system  in  time  At  is  split  into  two  stages. 

The  first  stage  models  the  change  of  internal  state  of  subsystems 
enclosed  in  the  cells,  assuming  that  the  particles  are  fixed;  collisions  of 
particles  (with  their  counterparts  in  a cell)  in  subsystems  {Cp...,^}  are 
simulated  independently  in  each  cell,  thus  the  particles  acquire  new  velo- 
cities. The  vector  c = {Cp...,^}  is  regarded  as  a state  of  Rats'  model. 

Let  4>(t,c)  be  the  density  of  the  probabilistic  distribution  of  the  state 
c(t);  then  the  governing  equation  of  this  model  ("Rats'  Master  Equation"  [24]) 
has  the  form 

' * 1<(U  ■ '22> 

where  R - Rats'  operator  of  collisions;  g.  = |c  -c  I , c.  and  c denote  the 
velocities  of  the  z-th  and  m-th  particles  upon  their  collision;  doJm  - a 
differential  section  of  elastic  dissipation  of  a pair  of  particles  (c? ,cm) ; 
a normalizing  parameter  V is  determined  by  the  choice  of  measurement  units 
and  it  can  be  interpreted  as  a cell  volume. 
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If  we  introduce  the  distribution  functions 


fsCt,clt...,c  ) = 

s (N-s)  IVs 


N 

f<P (t,c)  n dc.  , 
i=s+l  1 


then  by  integrating  (22)  it  is  not  difficult  to  obtain 

3t  = ! ^2^t,ci,cP  " f2^t,Cl’C2^g12*^CT12*^C2 


which  coincides  with  the  Boltzmann  equation  having  a zero  convective  deri- 
vative when  satisfying  equality(21) . 

The  algorithm  for  the  realization  of  the  first  calculational  stage  of 
evolution  of  a space -inhomogeneous  model  corresponds  to  the  Monte  Carlo 
method  of  numerical  solution  of  Rats'  basic  equation  (22)  which  (unlike 
the  Boltzmann  equation)  is  linear. 

The  second  stage  models  a collisionless  transfer  of  particles  from  a 
particular  cell  to  any  neighboring  one  without  changing  the  internal  state 
of  the  subsystems ; their  interactions  with  the  control  volume  boundary  and 
body  surface  are  considered  as  well.  This  stage  corresponds  to  the  Monte 
Carlo  method  of  numerical  solution  of  the  Boltzmann  free -molecular  equation 
in  the  following  form 


3f 

at 


+ 


c L f 


0. 


(23) 


where  L is  a finite-difference  operator  approximating  the  derivative  a/3x; 
its  introduction  is  closely  related  to  an  incomplete  description  of  the 
system  state  in  coordinate  space. 

The  simplest  numerical  algorithms  of  the  method  described  [21-23] 
correspond  to  the  solution  of  time  explicit,  conventionally- stable  finite- 
-difference  schemes  of  first-order  accuracy,  respectively,  for  Rats' 
equations  and  the  Boltzmann  free -molecular  equation.  Herein  the  equation 
of  evolution  of  a model  gas  (Rft^) , c(t  )}  within  the  accuracy  of  satisfying 
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molecular  chaos  has  the  form  (a  one-dimensional  flow) 


A f A f 

a + # a 

At  x ’ Ax 


J(ffj)  - AtCx 


AQJ(ffi) 

AX 


(24) 


where  A^  / At  and  <1^/Ax  are  first-order  finite-difference  operators  approxi- 
mating derivatives  3/3t  and  3/3x,  respectively;  JCff-^)  designates  the 
right-hand  side  of  the  Boltzmann  equation.^  The  finite -difference  scheme 
given  is  conventionally- stable  and  it  approximates  the  Boltzmann  equation 
within  the  accuracy  of  0(At)  and  0(ax).  As  mentioned  above,  this  calcula- 
tion employs  incomplete  information  concerning  the  space  position  of 
particles . 

This  calculational  model  can  naturally  be  extended  to  the  cases  of  a 
two-  and  three-dimensional  space;  it  consists  of  a sequence  of  one-dimen- 
sional displacements  along  each  coordinate  axes.  This  corresponds  to  the 
splitting  of  a multidimensional  transfer  equation 


into  a sequence  of  one -dimensional  finite -difference  schemes. 

The  Boltzmann  equation  is  known  to  imply  a molecular  chaos  or  a statical 
independence  of  particles.*  The  same  premises  are  inherent  in  our  model  as 
in  the  Boltzmann  equation  but  without  the  assumption  of  molecular  chaos 
(statistial  independence).  In  our  model  statistical  dependence  of  the  par- 
ticles violates  molecular  chaos.  It  should  be  noted  that  the  inherent 
statial  independence  rests  upon  theoretical  and  physical  premises  and  does 
not  depend  upon  the  mesh  dimension  (it  exists  as  Ax  -*•  0 as  well) . 

* The  molecular  chaos  hypothesis  implies  that  particle  velocities  are  sta- 
tistically independent.  (M.N.  Kogan,  "Rarefied  Gas  Dynamics,"  M. , Nauka,  1967). 
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The  numerical  results  of  rarefied  gas  flow  reveal  that: 

1)  results  of  calculations  with  various  numbers  of  particles  in  a cell 
(e.g.  with  N = 3 and  N = 20)  practically  coincide; 

2)  these  results  are  in  good  agreement  with  the  solution  of  the  Boltz- 
mann equation  (Cheremisin's  and  Rykov's  data).  Therefore,  the  violation  of 
molecular  chaos  in  the  problems  involved  is  small  (though  statical  depen- 
dence exists,  it  is  weakly  manifested  in  rarefied  gas  problems,  and  appar- 
ently, it  can  be  neglected  here) . 

For  the  study  of  turbulence,  statistical  independence  is  of  crucial  sig- 
nificance and  we  expect  that  it  will  be  manifested  in  this  method  when  applied 
to  turbulent  flows. 

5.2  Simulated  Structure  of  a Shock  Wave 

The  model  was  tested  for  the  solution  of  a problem  dealing  with  the 
structure  of  a normal  shock  wave  in  a gas  consisting  of  elastic  spheres  in 
the  range  of  Mach  numbers  =1.25  to  4. 

Graphs  of  density  n(x) , longitudinal  temperature  T11  (x) , transverse 
temperature  T^(x)  and  static  temperature  T(x)  are  shown  in  Figs.  17  and  18 
for  M = 2 and  3.  The  unit  of  length  is  the  mean  free  path  of  molecules 
in  the  free  stream  flow.  The  value  of  At/Ax  was  chosen  to  insure  stability 
of  the  calculations.  The  average  number  of  particles  in  cells 
corresponding  to  the  oncoming  flow  is  NQ  = 15  to  20  (M  = 2)  and  NQ  = 12 
(M  = 3).  A comparison  is  made  in  the  figures  with  the  density  h(x)  and 
temperature  T(x)  obtained  by  direct  numerical  integration  of  the  Boltzmann 
equation  [35,36]  on  a network  Ax  similar  to  the  one  used  in  our  calculations 
(Ax  = 0.2  to  0.3) . 

Finally,  Figs.  19  and  20  show  that  the  results  of  the  calculations  are 
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6.  Numerical  Investigation  of  Some  Gas  Dynamics  Problems  by 
Net-Characteristic  Methods. 

The  manifold  problems  presently  confronting  computational  mechanics 
are  increasing  in  complexity,  which  in  turn  requires  the  improvement  of 
earlier  numerical  methods  and  the  creation  of  new  ones.  The  properties 
that  a numerical  method  is  to  be  endowed  with,  from  the  view-point  of  modem 
developments,  are  so  diverse  that  they  are  difficult  to  fully  implement  in 
one  single  method.  In  view  of  this,  complexes  of  numerical  methods  based 
upon  a unified  approach  should  be  available. 

A large  and  important  class  of  problems  is  described  by  multi -dimen- 
sional systems  of  equations  of  hyperbolic  type.  A fundamental  concept 
underlying  the  construction  of  numerical  approximations  for  such  systems 
requires  that  their  characteristic  properties  be  taken  into  account  in 
seme  form,  i.e.,  it  is  essential  to  impart  relevant  features  to  the  numeri- 
cal methods  of  their  solution  as  well.  Finally,  it  is  desirable  to  consi- 
der homogeneous  schemes  that  enable  calculations  through  discontinuities 
that  may  arise  in  the  evolution  of  the  solution,  that  allow  for  explicitly 
singling  our  some  (principal)  of  the  features  and  that  adequately  resolve 
their  boundary  conditions. 

6.1  Investigation  of  Difference  Schemes  for  a Model  Equation 

Positive  type  difference  schemes,  first  introduced  in  [37],  play  an 
important  role  in  the  solution  of  equations  of  hyperbolic  type.  Using  the 
method  of  indefinite  coefficients  and  the  characteristic  properties  of  equa- 
tions of  hyperbolic  type,  one  can  write  a rather  large  class  of  schemes  of 
this  kind  and  then  perform  a comparative  analysis,  in  particular,  with 
respect  to  the  magnitude  of  "approximation  viscosity"  (see,  e.g.,  [38]). 
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These  schemes  possess  a number  of  indisputable  advantages,  the  main  ones 
are  the  absence  of  oscillations  in  the  numerical  simulation  of  non-smooth 
solutions,  and  the  construction  of  efficient  algorithms  for  the  calculation  of 
boundary  points.  However , since  they  are  all  of  first-order  accuracy,  the 
utilization  of  even  the  best  of  the  schemes  [39,40]  (having  a minimal 
"approximation  viscosity")  requires  in  some  cases  a great  number  of  grid 
points  in  the  difference  net  and,  consequently,  a voluminous  number  of 
calculations.  Some  schemes  of  higher  order  accuracy  can  be  constructed  by 
a similar  method  in  which  difference  approximations  are  treated  as  compo- 
nents of  a linear  space  of  indefinite  coefficients. 

The  basic  idea  of  the  method  is  illustrated  by  the  simple  wave  equation 
v + xvx  = 0,  x = const.,  X > 0 (25) 

discretized  in  a pattern  typical  of  explicit  schemes  comprising  six  points 
(Fig.  21) 


(i  » » ^n-2^  * * * * » ^n*  xm+2^  » 


(26) 


Linear  difference  schemes  written  for  the  pattern  (26)  in  the  following  form 


v£+1  = Z av  v£*v  = Z a vn 
m n,v  w p p=-2  w m+lJ 


(27) 


when  substituted  into  (25)  yield 

Z (p  - v o)  aV  = - a 
p,v  11 


(28) 


Z a =l,o=  X-r/h  > 0 , 
P,v  v 


where  coefficients  a j and  are  excluded. 

Any  point  in  the  space  of  the  coefficients  that  remain  indeterminate 
a = [a^,  a0»a2*  22a, b)  gives  rise  to  a difference  scheme  of  first- 
-order  accuracy  for  eq.  (25);  when  this  point  is  inclosed  by  polyhedrons 
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AjA3 A6  with  0 < a < 1 (fig.  22a)  or  A^A^Ag  with  1 < a < 2 (fig.  22b) 

it  gives  rise  to  difference  schemes  of  positive  approximation  for  which  all 
coefficients 

av  > 0 . (29) 

P - 

Points  A^  with  0 < a < l (scheme  [39])  and  A2  with  1 < a < 2 correspond  to 
difference  schemes  with  the  least  "approximation  viscosity",  i.e.,  with  the 
smallest  value  of  the  coefficient  of  v in  the  first  differential  approxi- 

AA 

mation  of  (27)  [38]. 

Plane  B1,...,Bg  (fig.  22) 

Z (ii-va)V  = or  aQ  = 3(a_2  + ^ + 1‘°2  (30) 

P,v  p 

for  the  pattern  (26)  constitutes  a two-parameter  family  of  difference 
schemes  of  second-order  accuracy  for  the  solutions  of  (25).  In  figs.  22a, b 
the  boundaries  of  the  region  of  stable  schemes  with  an  approximation  order 
higher  than  the  first  are  shown  by  ticked  lines  on  plane  (30) . It  is  seen 
that  with  0 < a < 2 this  region  is  not  empty.  It  is  shown  in  [41]  that  there 
are  no  difference  schemes  of  the  form  (27)  having  a second  order  approxima- 
tion (a  higher  order  accuracy  as  well,  see  [38])  which  satisfies  the  con- 
straints of  eq.  (29),  i.e.,  plane  (30)  does  not  intersect  a closed  polyhe- 
dron (28),  (29). 

A straight  line  being  an  intersection  line  of  plane  (30) 

and  the  plane 

Z (a-va)3  av  = -a3  , (31) 

p v p 

incorporates  a one -parameter  family  of  schemes  of  third-order  accuracy. 

With  0 < o < 1 the  segment  (fig.  22a)  contains  stable  schemes  of 
third-order  accuracy.  Point  on  the  pattern  given  is  the  only  difference 
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scheme  of  fourth-order  accuracy.  The  familiar  Lax-Wendroff  scheme  is  indi- 
cated by  point  B,-  in  figs.  22a, b. 

6.2  Construction  of  Positive-Type  Difference  Approximations. 

Evidently,  various  kinds  of  oscillations  of  nonsmooth  solutions  that  ^ 

are  observed  in  familiar  difference  schemes  with  an  approximation  order  | 

higher  than  the  first,  and  that  are  not  present  in  schemes  of  positive  type 
approximation,  are  due  to  the  fact  that  some  of  the  coefficients  in  j 

difference  expressions  of  the  form  (27)  are  negative  for  schemes  of  i 

higher  order  accuracy.  It  is  natural  to  suppose  that  the  behavior  of  a 
particular  scheme  for  nonsmooth  solutions  (the  amplitude  and  character 

j 

of  oscillations)  is  determined  by  the  distance  of  a point  (corresponding  to 
this  scheme  in  the  space  of  indeterminate  coefficients)  from  the 
region  of  difference  schemes  with  positive  approximation  (from  polyhedrons 
AjA3...A6  (fig.  22a),  A2A^A4A6  (fig.  22b)  etc.).  From  this  observation, 
it  is  proposed  that  the  "non-monotonicity"  of  difference  schemes  be  charac- 
terized by  the  value 


where  a(a^}  is  the  set  of  coefficients  in  eq.  (27)  corresponding  to  the 
difference  scheme  in  question  (of  approximation  order  higher  than  the  first 
situated  on  plane  (30)),  and  is  the  set  of  coefficients  in  eq.  (27) 
corresponding  to  the  vertex  of  the  polyhedren  prescribed  by  eqs.  (28),  (29) 
which  describe  difference  schemes  of  first-order  accuracy  of  positive  type 
approximation  (point  Aj  in  fig.  22a  with  0 < a < 1 and  point  A2  in  fig.  22 
with  1 < o < 2 etc.).  That  is,  schemes  of  first-order  accuracy  should  be 
constructed  on  three-point  patterns  incorporating  the  point  (tn+\  x^)  that 
is  to  be  calculated  and  two  nodes  of  the  difference  molecule  along  tn  which 
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just  brackets  the  characteristics  dx  = xdt  (e.g.,  xm_2>xm-i  ^8*  21).  •.  I 

For  difference  schemes  for  the  solution  of  (25)  of  second-order  ac- 
curacy, the  coefficients  are  found  by  a conventional  geometric  construc- 
tion in  the  a-space  (a  = {a^ag.c^})  of  the  point  of  intersection  with  i 

plane  (30)  of  the  normal  drawn  from  point  A.  With  0 < cr  < 1 (a^  = (0,1- a, 0}) 
this  procedure  results  in  point  in  fig.  22a,  a difference  scheme  of 
second-order  accuracy.  With  1 < a < 2 such  a scheme  corresponds  to  point 

/ 

Bg  in  fig.  22b.  (aA={o-l,  0,0}). 

The  difference  scheme  of  third-order  accuracy  with  the  smallest  value 
of  y which  is  stable  for  0 < a < 1 is  found  by  constructing  the  point  of 
intersection  of  the  normal  drawn  from  point  aA  = {0,l-a,0}  with  the 
straight  line  ((30),  (31)),  point  C3  in  Fig.  22a. 

The  calculation  of  the  simplest  modelling  problems  reveals  that  dif- 
ference schemes  of  second  and  third-order  accuracy  constructed  in  this 
manner  have  the  shortest  amplitude  of  oscillations  with  fast  damping  as 
compared  to  other  schemes. 

Bearing  in  mind  that  eq.  (25)  is,  in  essence,  an  ordinary  differential 


equation  along  the  characteristic  dx  = xdt 

dv_r>  d _ 3 3 

dt-0’  at"3t+X*37 

and,  consequently,  at  point  (tn+\xm)  the  exact  value  is 


n+i  f.m 

vm  ’ v(t  -x>  • x ’ =Sn  • 


so  that  the  difference  expressions  in  the  right-hand  part  of  (27)  are  es- 
sentially interpolation  formulas  for  calculating  v(tn,x) . From  this  view- 
point the  right-hand  part  of  eq.  (27)  for  schemes  possessing  the  smallest 
value  of  y may  be  treated  as  interpolation  polynomials  of  the  second  and 
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third-order  for  the  function  v(tn,x);  these  polynomials  differ  from  a 
piecewise -linear  interpolation  relevant  to  difference  schemes  of  first - 
-order  accuracy  with  positive  approximation.  A more  detailed  description 
of  the  approach  in  hand  including  "implicit"  difference  molecules  will  be 
given  in  a forthcoming  publication. 

6.3  Positive-Type  Difference  Schemes  for  Model  Equations 
The  construction  of  difference  schemes  of  sections  6.1  and  6.2  is  rea- 
lized by  a separate  compatibility  condition  along  characteristics  dx=X^dt 
of  a one-dimensional  system  of  equations  of  hyperbolic  type 

Ut  ♦ AUx  = f , A = a'1Afl  , (33) 

which,  in  a conventional  way,  is  reduced  to  canonical  form 


wiUt  + Xiwiux  = wif  > i = 1.....I  • (34) 

Here  U = {Uj,...,Uj}  is  the  vector  of  unknown  functions;  f = { fj , . . . , f j ) - 
is  a vector-column  of  right-hand  parts;  A = is  a diagonal  matrix  from 
eigenvalues  of  matrix  A;  52  = {w^}  is  a nonsingular  matrix  whose  lines  are 
linear -independent  eigenvectors  w^  of  the  matrix  A. 

If  matrix  A has  fixed  components  and  f = 0,  difference  schemes  (27) 
are  generalized,  in  an  obvious  manner,  for  the  case  of  system  (9) 
u"+1  = 2 DV  l rv  • DV  = 52_1  Av  52  , 


m 


u,v 


y m+y 


where  A^|  = { (a^) are  diagonal  matrices.  The  general  quasilinear  case  of 
system  (33)  (A  = A(t,x,v),  f = f(t,x,v))  requires  the  development  of  a pro- 
per method  of  integration  of  "ordinary"  differential  equations  (34)  that  takes 
into  account  the  dependence  of  X^,w^,f  upon  the  unknown  solution  v and  inde- 
pendent variables  t,  x.  Different  approaches  may  be  used  for  this  purpose. 

One  example,  intended  for  the  construction  of  explicit  difference  schemes 
using  the  pattern  of  type  (26) , is  the  Runge-Kutta  method  applied  in  [43] 
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on  the  pattern  (26)  a one-parameter  family  of  difference  schemes  of 
third-order  accuracy  was  constructed 

CS  ■ - -C  - * «tCi  ♦ 0/2  • (36> 

CB  • >fi  - 6'<C  - ♦ »*<♦!£  * C)/2  • (!7) 

(Fn+8_  1 

TIn+l  _ .n  __  3Tr  '•rm+l  m-l-*  , n+Bi  . t r vrm+l  ' m-lJ  . ,nn  . 


+ ♦“]  + 
mJ 


<-pn+B_  pn+g-j  / pH  pn  ■> 

if1*1  = if1  + m+i  m-.iJ  + *n+ei +i.[.  + ^ + 

m m TL  2h  *m  J TL  2h  V ^ 

+ t (ij  7-2fJ\,  +2F"  t-F11  -)/12h+  g(Un  ~-4Un  .1+6Un-4Uni,+u"  .)  , 
m-2  m+1  m-1  m-2"  m-2  m-1  m m+2  m+2J  ’ 

a = 1/3  , 6 = 2/3. 

It  was  also  shown  there  that  to  insure  stability  of  these  schemes  the 

value  of  the  parameter  g is  subject  to  the  condition: 

-1/8  < g < a*(a*-4)/24  , a*  = T max|x.|/h.  (39) 

m,i 

According  to  the  analysis  carried  out  in  sec.  6.1,2  for  linear  equations, 

it  is  preferable  to  choose  a matrix 

fi'1  • Gn  , G = (g.)  , n = (w.)  , 

1 (40) 

8j  = l0il  * (5lCTil  - 24)/152  , ai  = x^r/h,  i = 


instead  of  a scalar  factor  in  the  last  term  of  (38)  which  is  the  same  for 


all  characteristics.  In  equation  (39)  and  (40)  x^  are  the  eigenvalues  of  the 
matrix  A * 3F/3U,  a is  a matrix  derived  above  from  the  eigenvectors,  and  G 
is  a diagonal  matrix.  The  modification  just  cited  of  the  scheme  (36-38)  for 
a linear  case  is  the  scheme  discussed  in  section  6.2  corresponding  to  the 
smallest  value  of  y We  should  also  note  that  the  last  term  in  (38),  modified 
in  accordance  with  (40) , is  more  conveniently  written  in  divergence  form 
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(41) 


(fi_1Gfi)n  , • (-U"  n 3U11  .+U11  _)  - 
v -'m+?s  v m+1  m m+1  m+2' 

-(n'-W*  . • (-if  7+3lf  ,-3lf+uFJ  .) 

'm-h  m-2  m-1  m m+1 

On  the  same  pattern  of  (26)  there  can  be  generalized  for  the  case  of 

quasilinear  systems  (35) , difference  schemes  of  second-order  accuracy  which 

are  "closer"  to  schemes  with  positive  approximation,  for  example,  in  the  form 

of  a two-step  system  with  a predictor  - the  Lax  scheme  with  a = 0.5  in  (36) 

and  a corrector 

C1  ■ f-Ci  * K - 3Ci  * C2>  - 
- (“‘Vt  «C-2  * 3Ci  - K * CP  ♦ f42’ 

+ if  - to**  - ^b/h  + + *n+h/2 

m v m+h  m-V  v m+?s 

which  is  stable  providing 

= t max  | X - | /h  < 2.  (43) 

m,i 

Here  A = { (a  -) . } , A?  = {(a?).}  are  diagonal  matrices  with  elements 


(a2h  = 

3(1  - . 

a.) (2  - 

aJ/19,  ] 

1 for  1 < a.  < 2 

(»-2>i 

= (a_2) 

i-  a- 

■ o.) 

l - 

(a-2)i 

= 3(1  + 

a.)(2  + 

■ o.)/19,  1 

for  -2  < a.  < -1, 

Ca2)i  = 

(a.2)i 

- (1  + 

ai) 

1 

(a2)i  = 

(o-2)i 

II 

W 

Q 

(-1  - 1° 

^|/19  for  la.|  < 1. 

It  is  seen  that  satisfying  the  more  restrictive  condition  of  o*  < 1, 
rather  than  (43),  the  scheme  of  (36),  (42),  and  (44)  is  conservative.  The 
last  three  terms  in  equation  (42)  coincide  with  the  familiar  Lax-Wendroff 
scheme  [42]  and  can  be  replaced  by  other  modifications  of  analogous  schemes 
[45]. 
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As  an  example,  in  fig.  23a,  results  for  t/t  = 52  of  numerical  solutions, 

using  the  second-order  accurate  scheme  £36) , (42) , (44)  (scheme  II) , the 

Lax-Wendroff  scheme  [42],  and  the  MacGormak  scheme  [45],  are  compared  with 

the  exact  solution  (dashed  lines)  for  a problem  concerned  with  one -dimensional 

wave  motion  in  gas  resulting  from  initial  conditions:  v(0,x)  = 0,  p(0,x)  = 

= p(0,x)  = 2 for  x < 0;  p(0,x)  = p(0,x)  with  x > 0.  In  fig.  23b  a similar 

comparison  is  given  for  the  scheme  (36) -(38)  with  (a*  = 1)  choosing  g^  from 

2 2 

(40)  (scheme  III)  and  g^  = a*(o-*  - 4)/24  = -1/8.  In  the  same  figure  a scheme 
of  first-order  accuracy  with  positive  approximation  [46]  (scheme  I)  is  com- 
pared with  the  exact  solution  and  schemes  of  third-order  accuracy. 

In  this  problem,  as  well  as  in  the  linear  case,  one  observes  the  im- 
provement of  "oscillation"  properties  in  schemes  II,  III  as  compared  to 
other  schemes  of  second  and  third-order  accuracy.  Note  that  scheme  I as 
far  as  its  accuracy  is  concerned  is  quite  comparable  with  schemes  II,  III 
in  the  calculation  of  shock  waves;  however,  it  requires  a very  fine  difference 
net  for  calculations  through  contact  discontinuities. 

6.4  Applications  to  Gasdynamics  Problems 

Without  elaborating  here  upon  the  generalization  of  the  difference 
schemes  described  to  multidimensional  cases  it  should  only  be  noted  that, 
for  schemes  of  first-order  accuracy  with  positive  approximation  (27),  (28), 
the  construction  procedures  are  performed  rather  formally  (see,  e.g., 
[39,46,47]).  For  difference  schemes  of  higher-order  accuracy  there  also 
exist  various  efficient  approaches  described  in  detail  in  scientific  publi- 
cations. 

In  the  past  years  systematic  investigations  have  been  conducted  on  the 
aerodynamics  of  bodies  of  complex  geometry,  on  numerical  modelling  of  pro- 
blems in  plasma  physics,  and  on  dynamic  problems  in  the  theory  of  elasticity 
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utilizing  some  of  the  difference  schemes  discussed  above.  The  prominent 
features  in  such  problems  are  the  complex  structure  of  the  unknown  solution, 
the  presence  of  regions  with  large  gradients,  and  discontinuous  functions. 

Stable  schemes  of  first-order  accuracy  [47]  prove  to  be  sufficiently 
effective  for  the  numerical  solution  of  problems  that  involve  a com- 
paratively small  number  of  discontinuity  surfaces;  they  can  be  expli- 
citly defined  formulating  on  them  appropriate  boundary  conditions. 

As  an  example,  in  fig.  24a  are  shown  the  stationary  bow  shock  patterns 
in  the  plane  of  symmetry  of  a supersonic  three-dimensional  inviscid  flow  of 
a thermally  nonconducting  gas  (adiabatic  index  of  k = 1.4)  around  a spheri- 
cally blunted  cone  of  0^  = 10°  half-angle  which  has  a segmentally  capped 
base  (ec  = 35°).  The  overall  length  of  the  body  is  equal  to  5.5  radii  of  the 
spherical  nose  blunting.  The  Mach  number  of  the  incoming  flow  was  M=2,  with 
the  angle  of  attack  varied  from  0 to  180°.  Fig.  24b  presents  the  correspon- 
ding pressure  distributions  along  this  body.  Other  problems  of  three- 
-dimensional  supersonic  flow  of  a radiating  gas  around  blunt  bodies,  solved 
by  means  of  the  numerical  method  of  [47]  under  the  assumption  of  thermo- 
chemical equilibrium,  are  described  in  detail  in  [48-52]. 

Solutions  to  problems  containing  singularities  (discontinuities)  within 
the  integration  region  are  obtained  (with  acceptable  accuracy)  without 
singling  them  out  explicitly  by  the  introduction  of  conservative  elements 
(see,  e.g.,  [38,46])  in  the  difference  scheme  of  [47].  A number  of  one-  and 
two-dimensional  problems  involving  the  interaction  of  laser  light  with 
matter  are  considered  by  this  scheme.  As  a typical  example  of  this  kind 
(taken  from  [46])  fig.  25  presents  at  a time  instant  t * 10  10  sec.,  appro- 
priate to  the  termination  of  the  impulse,  the  isochores  p/pQ  * const,  (fig. 
25a)  and  isotherms  in  Kev  (the  electronic  temperature  Te  ■ const,  solid 


k 

I 

i 

| curves,  and  the  ionic  plasma  temperature  Ti  = const.  - dashed  curves,  fig. 

25b)  in  the  interaction  of  a synmetrical  impulse  of  laser  light  of  energy 
E = 300j  with  a spherical  envelope  of  variable  radius.  The  initial  distur- 
; bances  of  the  envelope  are  located  at  its  hali -radius  and  the  following 

physical  processes  are  taken  into  account;  the  absorption  of  outer  laser 
radiation,  nonlinear  electronic  thermal  conduction,  and  electron- ion  colli- 
sional  relaxation.  The  terms  in  the  energy  equations  for  electronic  and 
ion  components  related  to  thermal  conduction  and  energy  exchange  between  the 
components  are  approximated  in  an  implicit  way.  The  direct  application  of 
the  schemes  of  [47]  for  the  solution  of  similar  problems  resulted  in  a sig- 
; nificant  violation  of  the  integral  balance  of  mass,  momentum  and  energy,  so 

that  the  solution  in  the  vicinity  of  a nonstationary  shock  wave  moving  along 
a cold  background  is  quite  unsatisfactory. 

The  last  example  of  the  net  characteristic  methods  is  shown  in  Fig.  26, 
an  instantaneous  picture  of  the  distribution  of  components  a , a , a of 
a stress  tensor  in  an  elastic  layer  of  finite  thickness  at  the  instant 
t = t*  = 0.029.  They  are  produced  by  the  nonstationary  loading,  P,  on  part 
of  the  upper  boundary  of  the  elastic  layer  which  is  supported  on  a perfectly 
plane  rigid  base.  Boundary  conditions  in  this  problem  are  given  as  follows: 
on  the  portion  AB  of  the  upper  boundary, 

v (t,x,l)  = - -J-  /tdt[P+2/®a  (t,x,l)dx]  see  the  inset  of  Fig.  26c, 

y n ° A yy 

: 0 

axyft,x,1)  “ 0; 

on  the  other  portion  of  the  upper  boundary, 

(t,x,l)  - a (t,x,l)  * 0, 

77  ' 

an  the  lower  boundary,  y * .5, 

yt,x,0.5)  - axy(t,x,0.5)  = 0. 

AO  is  a plane  of  symmetry  and  m , the  mass  density  of  the  elastic  layer,  is 

I 
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a constant. 


In  conclusion,  the  author  would  like  to  express  his  sincere  gratitude 
to  V.V.  Demchenko  and  I.B.  Petrov  who  were  helpful  in  obtaining  the  numerical 
solutions  of  some  of  the  problems  described. 
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M„*  0.9 


Fig.  1 Lines  of  constant  Mach  number  (isotachs)  in  transonic  flow  past  a 

two-dimension  24  per  cent  circular  arc  profile;  critical  Mach 

number,  M * = 0.65. 
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Fig.  2 Isotachs  in  transonic  flow  past  a 24  per  cent  axisymmetric  body  (gen- 
erated by  revolution  of  a circular  arc  profile) ; critical  Mach  number. 


Supersonic  flow  patterns  around  a thick  circular  disk 


upstream. 


Fig.  6 Supersonic  flow  patterns  arourv'  a sphere. 
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300  J implosion  of  a pellet  (---  initial  inner  by  a nonstationary  loading  (shown  in  inset 

and  outer  boundaries  of  pellet).  of  an  elastic  layer  of  unit  thickness. 

Corresponding  isotherms  (Kev) ; electron 

temperature,  ---  ionic  plasma  temperature. 
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